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Research  Objectives 

The  primary  goal  of  this  work  was  to  understand  in  detail  the  relative  contributions  of  water 
column  optical  properties,  bottom  morphology,  bottom  material  reflectances,  bottom 
bidirectional  reflectance  distribution  functions  (BRDFs),  and  external  environmental 
conditions  on  remote-sensing  reflectances  in  optically  shallow  waters.  A  secondary  goal  of 
the  work  was  to  participate  in  various  collaborations  with  other  investigators  on  oceanic 
optics  problems  of  Navy  interest. 

Various  methodologies  are  now  under  development  for  the  extraction  of  environmental 
information  such  as  water-column  absorption  and  scattering  properties,  bottom  depth,  and 
bottom  type  from  remotely  sensed  hyperspectral  imagery  obtained  in  optically  shallow 
waters.  Regardless  of  the  methodology  used,  errors  in  measured  or  predicted  hyperspectral 
remote-sensing  reflectances  Rn(\)  will  degrade  our  ability  to  extract  information  from  the 
spectra.  The  potential  errors  therefore  must  be  understood. 

A  “look-up-table”  (LUT)  methodology  for  extraction  of  environmental  information  from 
measured  /?re  spectra  is  under  development  with  separate  funding  (Mobley  et  al.,  2004).  That 
technique  relies  on  matching  computed  and  measured  spectra.  To  assess  the  potential 
errors  in  the  spectrum  matching,  and  to  ascertain  where  additional  effort  should  be  expended 
in  improving  the  underlying  LUT  databases,  it  is  necessary  to  know  when  and  how  each 
potential  source  of  error  in  computed  spectra  comes  in  to  play. 

Other  research  performed  under  this  contract  include  investigations  of  radiometer  self 
shading,  underwater  visibility,  analytical  modeling  of  remote  sensing  reflectances  in  both 
deep  and  shallow  Case  2  waters,  semianalytical  methods  for  inversion  of  the  radiative 
transfer  equation,  and  the  use  of  neural  networks  and  multilinear  regression  for  inversion  of 
ocean  color  spectra. 
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Approach 

In  my  primary  research,  I  used  a  combination  of  HydroLight  (www.hydrolight.info;  Mobley 
and  Sundman,  200 1  a,  200 1  b)  and  Monte  Carlo  numerical  modeling  to  quantify  how  various 
sources  of  error  influence  predicted  Rn  spectra.  For  example,  one  can  expect  that  water- 
column  absorption  and  scattering  properties  (including  phase  function  effects)  will  be  less 
(more)  important  for  shallow  (deep)  waters,  and  that  bottom  properties  (BRDF,  material 
reflectance)  will  be  more  (less)  important  in  shallow  (deep)  waters.  However,  the  interplay 
of  these  error  sources  is  complex  and  simple  rules  for  error  analysis  are  hard  to  develop. 
Detailed  numerical  simulations  and  validation  with  observational  data  are  needed  for  full 
understanding. 

In  my  secondary  work,  I  used  both  HydroLight  and  my  Backward  Monte  Carlo  Three 
Dimensional  (BMC3D)  code  to  carry  out  the  needed  simulations. 

Work  Completed 

The  primary  work  quantifying  sensor,  sky,  water-column,  and  bottom  effects  on  retrieval  of 
environmental  information  from  remote-sensing  reflectances  was  presented  at  Ocean  Optics 
XVII,  held  in  Fremantle,  Australia  in  October  2004.  Those  results  are  now  being  prepared 
for  submission  as  a  refereed  journal  article. 

Additional  work  completed  and  published  under  this  contract  includes  the  following: 

•  A  study  of  self-shading  of  in-water  radiometers  (Leathers,  Downes,  and  Mobley, 
2004;  this  paper  is  attached  as  Appendix  A),  which  included  shallow  water  effects. 
This  work  used  the  BMC3D  Monte  Carlo  code  previously  developed  with  ONR 
funding  under  the  CoBOP  program. 

•  Development  of  an  analytical  model  for  prediction  of  remote-sensing  reflectances  in 
both  deep  and  shallow  waters  Case  2  waters  (Albert  and  Mobley,  2004;  this  paper  is 
attached  as  Appendix  B),  which  used  my  HydroLight  model. 

•  A  study  of  bioluminescence  and  underwater  visibility  (Johnsen,  Widder,  and  Mobley, 
2004;  this  paper  is  attached  as  Appendix  C),  which  used  a  modified  version  of  my 
BMC3D  Monte  Carlo  code  for  computation  of  the  needed  Point  Spread  Functions. 

•  A  comparison  of  neural  networks  and  regression  algorithms  for  inversion  of  remote¬ 
sensing  reflectances  (Dransfeld,  Tatnall,  Robinson,  and  Mobley;  this  paper  is 
attached  as  Appendix  D),  which  used  HydroLight  to  generate  test  data. 

•  Investigation  of  shape-factor  models  for  use  in  retrieving  inherent  optical  properties 
from  remote  sensing  reflectances  (Hoge,  Lyon,  Mobley,  and  Sundman,  2003;  this 
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•  Investigation  of  shape-factor  models  for  use  in  retrieving  inherent  optical  properties 
from  remote  sensing  reflectances  (Hoge,  Lyon,  Mobley,  and  Sundman,  2003;  this 
paper  is  attached  as  Appendix  E),  which  used  HydroLight. 

•  Completion  of  a  technical  report  documenting  the  Monte  Carlo  techniques  used  in 
several  of  my  studies  (Leathers,  Downes,  Davis,  and  Mobley,  2004;  this  report  is 
available  on  CD,  published  by  NRL-DC). 


Results 

Consider,  as  one  example  of  my  primary  investigations,  the  LUT  methodology  for  retrieval 
of  bathymetry  in  optically  shallow  waters.  The  HydroLight-generated  Rts  spectrum  database 
used  to  match  the  image  Rn  spectra  was  created  with  various  sets  of  inherent  optical 
properties  (IOPs,  namely  the  absorption  a,  scattering  b,  and  backscattering  bh  coefficients). 
Suppose  the  database  includes  RIS  spectra  corresponding  to  a  and  b  coefficients  that  are 
representative  of  the  imaged  water,  but  does  not  have  in  it  Rn  spectra  corresponding  to  the 
backscatter  fraction  B  =  b/bh  values  that  occurred  in  nature  at  the  time  the  image  was 
acquired.  The  database  Rts  spectra  will  then  be  a  mismatch  for  the  image  RIS  spectra,  and 
errors  in  the  retrieved  bathymetry  will  result.  We  thus  ask  how  important  it  is  that  the 
database  contains  IOPs  with  the  correct  B\  i.e.,  how  large  are  the  bathymetry  errors  if  the 
database  B  values  do  not  correspond  to  those  in  nature. 

Figure  1  shows  four  points  in  the  vicinity  of  Lee  Stocking  Island,  Bahamas,  for  which  water- 
column  IOPs,  bathymetry,  and  bottom  classification  are  known  from  field  observations  and 
are  well  described  in  the  existing  LUT  database.  These  points  correspond  to  shallow  (2  m 
depth)  and  deeper  (8  m)  sand  (highly  reflecting)  and  seagrass  (dark)  bottoms.  Figure  2  shows 
the  corresponding  PHILLS  spectra  and  the  closest-matching  spectra  from  the  LUT  database. 

Starting  with  the  “known  answers”  for  these  points,  I  generated  RIS  spectra  for  B  values 
ranging  from  0.005  to  0.05  (total  backscatter  fraction,  including  water  and  particles),  whereas 
the  correct  value  was  near  0.02.  The  depths  (2  m  and  8  m)  and  bottom  types  (sand  and  grass) 
were  held  the  same.  I  then  treated  the  newly  generated  RIS  spectra  as  though  they  were  image 
spectra  and  searched  the  database,  which  contained  only  spectra  with  B  near  0.02,  to  retrieve 
the  known  depths. 

Figure  3  shows  a  histogram  of  the  depths  retrieved  when  the  database  does  not  contain  the 
correct  water  column  B  values.  For  the  2  m  depth,  the  retrievals  always  gave  the  correct 
depth,  which  was  included  in  the  database  at  intervals  of  0.25  m.  In  other  words,  the 
mismatch  between  the  actual  (for  various  B  values)  and  database  spectra  was  not  enough  to 
trigger  selecting  a  bottom  spectrum  that  was  off  by  as  much  as  0.25  m.  For  such  shallow 
water,  the  water  column  backscatter  fraction  is  thus  unimportant.  For  the  8  m  depth,  the 
correct  depth  of  8  m  was  retrieved  80%  of  the  time  over  the  bright  sand  bottom,  but  there 
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was  a  spread  of  retrieved  depths  (7  to  9.5  m)  over  the  dark  grass  bottom.  This  shows  the 
importance  of  having  the  correct  backscatter  fraction,  all  else  being  the  same,  as  a  function 
of  bottom  depth  and  type. 

As  another  example  of  this  type  of  study,  I  generated  IOPs  with  various  a,  b,  and  bh  spectra: 
namely  for  5  sets  of  a-aw  5  sets  of  b-bw,  and  5  backscatter  fractions,  for  a  total  of  125 
combinations  of  IOPs.  Figure  4  shows  the  corresponding  Rn  spectra  for  the  absorption 
perturbations  (with  the  scattering  and  backscattering  spectra  held  constant  at  their  baseline 
values).  Figure  5  shows  the  histogram  of  depth  retrievals  for  the  125  combinations  of  IOP 
perturbations.  We  again  see  that  the  depth  retrievals  are  quite  good.  Even  for  the  dark 
seagrass  bottom  and  8  m,  for  which  the  IOPs  are  most  important,  the  depth  is  retrieved  to 
within  ±1  m  in  85%  of  the  cases. 

Similar  studies  were  carried  out  to  assess  the  importance  of  the  bottom  BRDF  (Lambertian 
vs  non-Lambertian  BRDFs),  bottom  irradiance  reflectance,  bottom  slope,  rippled  bottoms 
(as  opposed  to  smooth  flat  bottoms),  sky  conditions,  sensor  noise,  and  the  like.  Without 
going  into  the  details  of  each  aspect  of  the  study,  the  overall  conclusions  from  these 
investigations  can  be  summarized  as  follows: 

•  Retrievals  of  depth,  bottom  reflectance,  and  water-column  IOPs  are  not  degraded  by 
typical  amounts  of  sensor  noise  (e.g.,  random  noise,  spikes,  dropoff  near  400  nm). 

•  Systematic  offsets  in  measured  Rrs  do  degrade  retrievals,  but  such  offsets  are  often 
easy  to  identify  and  correct. 

•  Non-Lambertian  bottoms  (not  included  in  the  present  LUT  database)  can  cause 
depth-retrieval  errors  of  -10%  (i.e.,  std  dev  of  retrieved  depths  /  true  depth  -0.1). 

•  Sun  angle  (30  to  60  deg)  and  off-nadir  viewing  direction  (out  to  -30  deg)  are  not 
critical  for  LUT  retrievals. 

•  Random  noise  in  IOPs  does  not  degrade  retrievals. 

•  Systematic  perturbations  in  IOPs  can  cause  depth  errors  of  -10%  for  darker  and 
deeper  bottoms. 

I  am  currently  completing  the  study  of  the  effects  of  environmental  conditions  on  remote¬ 
sensing  reflectances  and  the  implications  for  LUT  retrieval  of  bathymetry  and  bottom 
classification.  This  work  was  described  in  a  talk  at  the  Ocean  Optics  XVII  conference.  A 
journal  article  presentation  of  the  full  suite  of  studies  is  in  preparation. 

The  results  of  the  secondary  studies  can  be  seen  in  the  attached  reprints  of  the  resulting 
journal  articles  and  therefore  will  not  be  discussed  here. 
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Impact  and  Applications 

Hyperspectral  imagery  is  increasingly  used  for  a  wide  range  of  problems  from  mapping  and 
monitoring  seagrass  beds  and  coral  reefs  to  remote  sensing  of  bathymetry  and  bottom 
classification  for  military  applications.  For  quantitative  analysis  of  hyperspectral  imagery 
it  is  necessary  to  have  calibrated,  accurate  hyperspectral  reflectance  spectra.  This  need  in 
turn  makes  it  necessary  to  evaluate  in  detail  the  various  sources  of  error  in  such  spectra.  By 
quantifying  various  error  sources,  this  work  provides  guidance  as  to  where  additional  effort 
should  be  expended  to  improve  measurements  and  models  used  in  the  analysis  of 
hyperspectral  data. 

The  secondary  investigations  contribute  to  our  overall  knowledge  of  underwater  visibility 
and  to  the  development  of  a  wide  range  of  tools  for  extraction  of  environmental  information 
from  remote-sensed  ocean  color  data  in  a  variety  of  waters  (deep  or  shallow,  Case  1  or  Case 
2).  All  of  these  topics  are  of  great  Navy  relevance,  as  well  as  of  interest  to  the  optical 
oceanography  research  community  in  general. 


Collaborations  and  Related  Projects 

This  work  used  data  sets,  imagery,  and  models  (namely  PHILLS  imagery  and  the  BMC3D 
computer  code)  previously  obtained  or  developed  during  the  ONR  CoBOP  program.  This 
work  directly  contributed  to  my  separately  funded  work  on  developing  the  look-up-table 
methodology  for  extraction  of  environmental  information  from  hyperspectral  imagery. 

My  primary  work  was  done  in  close  collaboration  with  Paul  Bissett,  Dave  Kohler,  Mubin 
Kadawala,  Bhavesh  Goswami,  and  others  at  the  Florida  Environmental  Research  Institute. 
PHILLS  imagery  and  other  collaborations  were  provided  by  Curtiss  Davis,  Robert  Leathers, 
Valerie  Downes,  Marcos  Montes  and  others  at  the  Naval  Research  Lab,  Washington  D.C. 

The  collaborations  on  the  secondary  work  can  be  seen  in  the  author  list  of  the  attached 
publications. 
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Figure  1 .  PHILLS  image  of  the  Adderly  Cut  area  near  Lee  Stocking  Island,  Bahamas.  The 
numbered  points  are  the  pixels  where  I  performed  detailed  calculations  to  evaluate  the 
various  sources  of  error  in  computed  remote-sensing  reflectances. 
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Fig.  2.  Measured  (PHILLS,  blue  curves)  and  closest-matching  LUT  database  (red  curves) 
spectra.  These  spectra  were  the  starting  points  for  the  study. 
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Figure  3.  Histograms  showing  the  errors  in  the  retrieved  bottom  depths  if  the  database  of 
remote-sensing  reflectances  does  not  contain  the  correct  water-column  backscatter  fraction. 
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Fig.  4.  The  Rrs  spectra  corresponding  to  perturbations  in  the  absorption  coefficient.  The 
upper  group  of  curves  is  for  sand  at  2  m  depth;  the  middle  group  is  sand  at  8  m;  the  bottom 
two  groups  are  seagrass  at  2  and  8  m. 
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Fig.  5.  Histogram  of  the  depth  retrievals  for  125  different  lOP  perturbations  (as  illustrated 
in  Fig.  4  for  the  absorption  perturbations.). 
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Abstract:  We  present  the  derivation  of  an  analytical  model  for  the  self¬ 
shading  error  of  an  oceanographic  upwelling  radiometer.  The  radiometer  is 
assumed  to  be  cylindrical  and  can  either  be  a  profiling  instrument  or  include 
a  wider  cylindrical  buoy  for  floating  at  the  sea  surface.  The  model  treats 
both  optically  shallow  and  optically  deep  water  conditions  and  can  be 
applied  any  distance  off  the  seafloor.  We  evaluate  the  model  by  comparing 
its  results  to  those  from  Monte  Carlo  simulations.  The  analytical  model 
performs  well  over  a  large  range  of  environmental  conditions  and  provides  a 
significant  improvement  to  previous  analytical  models.  The  model  is 
intended  for  investigators  who  need  to  apply  self-shading  corrections  to 
radiometer  data  but  who  do  not  have  die  ability  to  compute  shading 
corrections  with  Monte  Carlo  simulations.  The  model  also  can  provide 
guidance  for  instrument  design  and  cruise  planning. 
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1.  Introduction 

When  a  radiometer  is  positioned  to  measure  upwelling  light,  its  shadow  decreases  the 
magnitude  of  the  local  light  field,  causing  measured  values  of  upwelling  radiance,  and  hence 
remote  sensing  reflectance,  to  be  too  low.  The  magnitude  of  the  error  depends  on  wavelength, 
sensor  size,  water  turbidity,  and  illumination  conditions  [1],  and  in  shallow  water  it  also 
depends  on  the  water  depth  and  seafloor  reflectivity  [2].  The  self-shading  error  has  a  direct 
effect  on  the  measurement  of  atmospheric  optical  thickness  and  on  remote  sensing  sensor 
calibration.  Furthermore,  because  the  magnitude  of  the  shading  error  is  wavelength 
dependent,  algorithms  that  depend  on  the  spectral  shape  of  upwelling  radiance  to  determine 
water  properties,  water  depth,  or  water  characteristics  [3-7]  will  also  return  erroneous  results. 
Therefore,  proper  ocean  optics  protocols  dictate  that  shading  corrections  should  be  routinely 
applied  to  all  in  situ  upwelling  light  measurements  [8].  Unfortunately,  instrument 
manufacturers  do  not  typically  provide  self-shading  correction  algorithms  for  their  products. 

We  developed  a  Monte  Carlo  code  [2]  to  compute  the  self-shading  of  buoyed  radiometer 
data  collected  in  1999  and  2000  at  Lee  Stocking  Island  (LSI),  Bahamas  [5-7].  The  self- 
shading  of  upwelling  radiometers  has  also  been  investigated  by  other  researchers  [1,  9-11]; 
however,  none  of  these  investigations  consider  optically  shallow  waters  such  as  those  at  LSI 
nor  do  they  account  for  the  additional  shading  caused  by  a  flotation  buoy  such  as  that  on  the 
Hyper-TSRB. 

The  objective  of  this  paper  is  to  present  a  new  analytical  self-shading  correction  model  for 
buoyed  and  unbuoyed  upwelling  radiometers.  We  evaluate  the  model  with  numerical 
simulations  and  then  discuss  the  model’s  uses  and  limitations.  Although  Monte  Carlo 
computations  provide  a  more  accurate  method  for  making  self-shading  corrections,  we  believe 
that  analytical  or  semianalytical  models  are  more  likely  to  be  implemented  by  the  general 
community  because  of  the  ease  of  which  they  can  be  disseminated  and  applied. 

2.  Model  derivation 

We  present  separate  analytical  shading  corrections  for  a  sensor  optically  far  from  the  seafloor 
and  a  sensor  close  to  the  seafloor  and  then  combine  the  two  into  one  general  correction 
algorithm.  We  derive  the  shading  correction  for  a  radiometer  far  from  the  seafloor  by 
considering  the  idealized  model  shown  in  Fig.  1,  which  is  a  slightly  generalized  version  of  the 
model  used  by  Gordon  and  Ding  [1].  A  small  sensor  with  infinitesimally  small  field  of  view 
(FOV)  is  positioned  at  a  distance  zs  below  a  shading  disk  of  radius  r.  The  shading  disk  can 
represent  the  bottom  of  either  the  sensor  head  (if  zs  =  0)  or  a  buoy  (if  zs  >  0).  The  solar 
illumination  is  taken  to  be  collimated  and  the  amount  of  in-water  scattering  by  water  is 
assumed  to  be  small  enough  to  not  significantly  disturb  the  collimated  nature  of  downwelling 
light.  The  goal  is  to  determine  how  much  the  measured  upwelling  radiance  is  reduced  by  the 
shadow  that  falls  across  the  sensor’s  line  of  sight. 

The  depth  to  which  the  shadow  lies  across  the  sensor’s  line  of  sight  is  [1] 

z0  =  r/tan^  ,  (1) 

where  the  in-water  solar  zenith  angle  6L*  is  related  to  the  above-water  solar  zenith  angle  6b  by 

[12] 

Ofto  =  sin'1  (sin0o/l  .338) .  (2) 
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Fig.  I .  Shading  of  an  upwelling  radiance  sensor  by  a  horizontal  disk  of  radius  r. 

The  depth  dependence  of  the  true  (i.e.,  unshaded)  radiance  Z,u*  can  be  expressed  as 

V  (z)  =  A*  (z,)exp[-*u  (*“*.)]»  (3) 

where  Ku  is  a  depth-averaged  value  of  the  diffuse  attenuation  coefficient  for  (z) .  Because 

a  negligible  amount  of  light  is  scattered  into  the  detector  from  within  the  shaded  region,  the 
measured  radiance  Lum  is  equal  to  the  radiance  at  depth  z0  attenuated  back  up  to  the  sensor 
depth  zs, 

Ium  =  Lu,(z0)exp[-a(z0-z,)].  (4) 

Using  Eq.  (3)  to  obtain  (z0)  and  substituting  this  into  Eq.  (4),  we  obtain 

V  =  V(z*)<*p[-(*u+«)(zo  "*.)]•  (5) 

For  collimated  light,  Ku  =  a/ COS6U;  substituting  this  and  Eq.  (1)  into  Eq.  (5)  gives 

Ium  =  (zs)exp[-&  a(r-z$  tan^  )],  tan^,  <  r/zs ,  (6) 

where 

k  =  (l/tan#^  +  1/sin^).  (7) 

The  fractional  shading  error  in  deep  water  is 

£  t  A,"  [l-expf-*  a(r-z, tanfl^)],  tan^  <r/z,| 

4'(z.)  {  0,  tan5>0w>r/zl  J 

For  self-shading  of  a  sensor  head  only  (i.e.,  zs  =  0),  Eq.  (8)  reduces  to 

Sy,=  [l-exp(-*ar)].  (9) 

Equation  (9)  is  the  same  expression  derived  by  Gordon  and  Ding  [1]  except  that  they  derived 
it  using  the  approximation  Ku  =  a  (instead  of  Ku  =  a/cos 6tw),  which  leads  to  k  =  2/tan Oo*.  The 
difference  between  k  =  2/tan6bw  and  k  =  (l/tan6bw  +  l/sin^)  is  only  significant  for  large 
values  of  6fcw. 

In  our  low-scattering  model,  the  shading  due  to  the  sensor  head  and  that  due  to  the  buoy 
are  never  additive.  At  small  solar  zenith  angles,  the  depth  of  the  buoy’s  shadow  below  the 
sensor  is  deeper  than  the  depth  of  the  sensor-head’s  shadow,  and  the  sensor  head  contributes 
no  shading  beyond  that  already  caused  by  the  buoy.  Conversely,  the  presence  of  the  buoy 
contributes  no  additional  shading  to  that  caused  by  the  sensor  head  at  large  solar  zenith 
angles.  The  transition  is  at  tan6fcw  =  (r8  -  rb)/zf,  where  r%  and  rb  are  the  radii  of  the  buoy  and 
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sensor  head,  respectively,  and  zg  is  the  distance  between  the  bottoms  of  the  two.  The  total 
shading  of  a  buoyed  radiometer  is  then  simply  the  larger  of  the  two  effects, 

=  jl-exp[-*a(rb-zstan0ow)],  tan^  <(rs -rb)/zs| 

l-exp[-*ars],  tan 6^  >(rs  -rb)/zs  j 


(10) 


Although  Fig.  1  shows  the  shading  disk  at  the  sea  surface,  Eq.  (10)  is  valid  for  submerged 
instruments  if  the  illumination  is  handled  correctly;  this  is  discussed  below.  Also  note  that  the 
derivation  of  Eq.  (10)  assumes  that  the  total  water  depth  Zbot  is  greater  than  the  depth  z0; 
otherwise,  =  1. 

An  expression  for  self-shading  close  to  the  seafloor  (i.e.,  when  the  upwelling  light  is  due 
primarily  to  bottom  reflection)  can  be  developed  using  the  model  shown  in  Fig.  2.  A  sensor 
with  finite  FOV  is  located  at  depth  zs  and  a  shading  disk  of  radius  rd  is  located  at  depth  zd.  (For 
the  case  of  self-shading  of  the  sensor  head,  zd  =  zg.)  The  total  water  depth  is  Zb*.  If  the  seafloor 
is  horizontal  then  both  the  sensor  FOV  and  the  shading  disk  project  circles  onto  the  seafloor. 
For  a  Lambertian  seafloor  and  a  sensor  that  responds  equally  to  light  from  all  directions 
within  its  FOV,  the  percent  shading  due  to  the  shadow  on  the  seafloor  approximately  equals 
the  fraction  of  the  FOV  circle  that  is  overlapped  by  the  shading  circle.  One  of  four  situations 
exists:  the  two  circles  do  not  overlap;  the  FOV  circle  is  completely  inside  the  shade  circle;  the 
shade  circle  is  completely  inside  the  FOV  circle;  or  the  two  circles  partially  overlap.  Let  the  x- 
axis  lie  in  the  horizontal  plane  such  that  the  positive-x  direction  is  looking  away  from  the  sun. 
The  shadow  on  the  seafloor  is  of  radius  rd  and  is  centered  at  xd  =  (tan#>H,)(zbot  -  zd).  The  FOV 
on  the  seafloor  is  centered  at  x  =  0  and  has  radius  rfov  =  (tan#ov)(>bot  ~  zs),  where  #ov  is  the 
FOV  half-angle.  The  fractional  shading  can  be  expressed  mathematically  as 


= 


0, 

1, 

2 


nrf . 


V*-r4 


!•{ 

**  J 

(*d 

(-'fev)>(*d-'i) 

(*d+'i)<rfcv 

(-'fo.)<(*d-'i)<'-f« 

and'ib.  <(*d+'i) 


,(H) 


where 

xd  —  tan#**  (zbot —  zd )  ,  r f0V  —  [tan#bv  (^bot —  £»)]>  <md  xmt 
The  integrals  in  Eq.  (1 1)  can  be  evaluated  analytically: 
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For  a  buoyed  radiometer,  we  will  approximate  the  overall  shading  error  close  to  the  seafloor 
with 

Sb  =  maxte  f  (12) 

where  ^  is  obtained  from  Eq.  (11)  with  rd  and  zd  taken  to  be  the  radius  and  depth  of  the  sensor 
head  and  3,  is  obtained  from  Eq.  (11)  with  rd  and  zd  taken  to  be  the  radius  and  depth  of  the 
buoy. 


#5125  -  515.00  US 
(C)  2004  OSA 


Received  25  August  2004;  revised  15  September  2004;  accepted  16  September  2004 
4  October  2004  /  Vol.  12,  No.  20  /  OPTICS  EXPRESS  4712 


Fig.  2.  Model  for  computing  the  amount  of  the  overlap  between  the  shadow  on  the  seafloor  and 
the  sensor  FOV. 


Now  that  we  have  developed  a  shading  model  for  a  radiometer  far  from  the  seafloor  [Eq. 
(10)]  and  a  shading  model  for  a  radiometer  very  close  to  the  seafloor  [Eqs.  (11)  and  (12)],  the 
final  step  is  to  construct  a  transition  from  one  to  the  other.  We  take  the  overall  fractional  error 
e  to  be  the  weighted  sum  of  the  two  models  based  on  the  relative  importance  of  the  signal 
from  the  seafloor  at  the  depth  of  the  radiometer,  i.e. 


(13) 


where  £*,  is  obtained  with  Eq.  (10),  €q  is  obtained  with  Eq.  (12),  Luw  is  the  portion  of  Z,u*  that 
originates  from  scattering  within  the  water  column  and  Z,„b  is  the  portion  of  that  originates 
at  the  seafloor  (Lub  +  £Uw  =  Lux).  The  ratios  (LUJ  Lu{)  and  (Lub  /  Lux)  can  be  computed  with  a 
numerical  radiative  transfer  code  such  as  Hydrolight  (Sequoia  Scientific,  Bellevue, 
Washington).  Alternatively,  (LUJ Lu*)  and  (Lub  / L»)  can  be  approximated  with 


4w(2.) 

4'(z.)  +(^>  aMo™  X ~K)cxV[~a X zt  \  ’ 


(14) 


where  zd  =  -  z8,  z=  (1+1/m>w),  and  Rb  is  the  seafloor  albedo  (0  <  Rb  <  1).  Our  derivation  of 

Eq.  (14)  has  been  omitted  for  the  sake  of  brevity.  The  value  of  Rb  can  be  chosen  either  by 
selecting  the  appropriate  value  from  published  bottom  reflectance  spectra  [6]  or  by  esimating 
its  value  from  a  model  of  the  measured  light  field  [13]. 

In  summary,  the  total  error  for  a  given  incident  illumination  direction  is  given  by  Eq.  (13), 
where  is  given  by  Eq.  (10),  sB  is  given  by  Eqs.  (1 1)  and  (12),  and  the  ratios  LUJLU{  and 
Iub/Lu1  are  either  computed  with  a  radiative  transfer  code  or  approximated  with  Eq.  (14). 

Equations  (10)  and  (11)  were  derived  for  collimated  illumination  and  low  in-water 
scattering.  However,  these  equations  can  be  applied  to  more  general  and  realistic  conditions 
by  computing  a  weighted  average  of  the  shading  errors  due  to  light  incident  on  the  instrument 
from  various  directions,  i.e.. 


e  =  ^wlsl  , 


(15) 
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where  the  weight  w*  of  each  angular  bin  is  proportional  to  the  power  of  light  incident  at  that 
directional  bin.  As  a  special  case  of  this  approach,  Gordon  and  Ding  [  1  ]  proposed  taking  the 
total  shading  error  to  be 

&skyf^  0  ~  f)  ^sun  >  0^) 

where  is  the  shading  error  computed  for  direct  sunlight,  is  the  shading  error  of  diffuse 
skylight,  and /is  the  fraction  of  the  total  downwelling  irradiance  that  is  skylight.  Numerically, 
the  diffuse  skylight  error  is  approximately  equal  to  that  for  direct  sunlight  at  Oq  =  35°  [1]. 
Alternatively,  values  of / can  be  computed  with  a  radiative  transfer  code,  such  as  Hydrolight. 
Example  values  of  /  for  use  in  Eq.  (16)  are  provided  by  Baker  and  Smith  [14].  For  a 
submerged  instrument  it  may  also  be  necessary  to  account  for  in-water  scattering  when 
estimating  the  angular  distribution  of  light  incident  on  the  instrument.  Furthermore,  in 
optically  shallow  waters  it  might  be  necessary  to  account  for  internal  reflection  at  the  sea 
surface.  In  the  last  two  cases,  one  should  estimate  6^  directly  for  use  in  Eqs.  (10)  and  (1 1) 
rather  than  with  Eq.  (2). 

3.  Assessment 

Gordon  and  Ding  [1]  used  Monte  Carlo  simulations  to  develop  a  semi-analytical  model  for  the 
self-shading  of  a  cylindrical  radiometer  far  from  the  seafloor.  This  model,  which  we  will 
hereafter  refer  to  as  the  GD  model,  is  simply  Eq.  (9)  with  the  value  of  k  replaced  with  a  value 
determined  empirically  from  their  Monte  Carlo  results.  This  value  of  k  depends  on  the 
illumination  conditions  and  on  the  size  and  angular  response  of  the  sensor.  The  GD  model 
was  verified  in  field  experiments  by  Zibordi  and  Ferrari  [9]  and  by  Aas  and  Korsbo  [10]  and 
is  endorsed  by  the  ocean  optics  protocols  for  SeaWiFS  validation  [8].  We  therefore  use  this  as 
a  starting  point  for  evaluation  of  our  more  general  analytical  model.  Table  1  shows  values  of  k 
obtained  three  different  ways:  with  k  =  2/tan 6U  (as  derived  analytically  by  Gordon  and  Ding 
[1]),  with  our  analytical  result  of  k  =  (l/tan6bw  +  l/sin^),  and  from  the  empirically  derived 
values  of  the  GD  model  for  a  point  sensor.  The  values  k  =  (l/tan^  +  l/sin^)  agree  with  the 
numerical  results  better  than  do  k  =  2/tan  6^,  especially  at  large  solar  zenith  angles.  It  seems 
likely  that  the  difference  between  the  GD  values  of  k  for  a  point  source  and  k  =  (l/tan6bw  + 
l/sin6tw)  is  due  primarily  to  the  uncertainty  in  their  Monte  Carlo  results.  This  assertion  is 
based  on  the  fact  that  the  numerical  values  are  slightly  larger  than  the  analytical  values, 
whereas  we  would  expect  the  opposite  because  the  analytical  model  ignores  the  effects  of 
scattering.  In  any  case,  our  analytical  values  of  k  lie  in  between  those  of  the  GD  model  for  a 
point  sensor  and  those  for  a  finite  sensor  (not  shown  in  Table  1)  and  therefore  provide  a  good 
general-purpose  model  for  the  self-shading  of  a  unbuoyed  radiometer  far  from  the  seafloor. 


Table  1 .  Values  of  k  [Eq.  (9)]  for  a  radiance  point  sensor 


ft (deg) 

k  =  2/tan  6U, 

£  =  (l/tan^  +  l/sin6bw) 

k  from  numerical 
simulations  [1] 

10 

15.28 

15.35 

16.58 

20 

7.56 

7.69 

8.43 

30 

4.96 

5.16 

5.54 

40 

3.65 

3.91 

4.18 

50 

2.86 

3.18 

3.39 

60 

2.36 

2.72 

2.84 

70 

2.03 

2.44 

2.48 

Our  assessment  of  the  rest  of  our  analytical  model  was  performed  with  simulations  from 
the  backward  Monte  Carlo  program  described  in  Leathers  et  al.  [2].  Many  of  our  numerical 
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computations  were  done  with  the  Hyper-TSRB  in  mind.  The  Hyper-TSRB  measures  above¬ 
water  downwelling  irradiance  (£d)  and  underwater  upwelling  radiance  (Lu)  at  many  spectral 
channels  from  400-800  nm,  from  which  remote  sensing  reflectance  spectra  can  be  computed. 
The  upwelling  radiance  sensor  has  a  radius  of  0.045  m  and  is  suspended  0.66  m  below  the 
sea-surface.  The  buoy  has  a  diameter  of  0.15  m,  and  the  bottom  of  the  buoy  sits  0.12  m  below 
the  sea  surface  and  0.54  m  above  the  upwelling  radiance  sensor.  The  angular  response  of  the 
radiance  sensor  was  provided  by  Satlantic;  the  effective  FOV  is  approximately  8.5  deg. 

Shown  in  Figs.  3  and  4  is  the  self-shading  error  for  the  Hyper-TSRB  in  optically  deep 
water  predicted  by  Eq.  (10)  and  for  the  sensor  head  only  predicted  by  Eq.  (9).  Also  shown  in 
Figs.  3  and  4  are  values  of  Hyper-TSRB  self-shading  error  computed  with  Monte  Carlo 
simulations  for  bla =  1 .  It  can  be  seen  that  the  buoy  adds  a  significant  amount  of  shading  for 
solar  zenith  angles  less  than  15°  but  contributes  no  additional  shading  beyond  that  of  the 
sensor  head  for  ft  >  15°.  Our  analytical  model  gives  excellent  agreement  with  the  numerical 
results  in  all  cases  except  for  the  combination  of  small  solar  zenith  angle  and  high  water 
scattering.  (Note  that  in  Figs.  3  and  4  the  value  of  b  is  large  when  a  is  large.)  This  gives  us 
great  confidence  in  our  analytical  model  for  most  practical  purposes.  It  should  be  noted  that 
the  case  of  ft  <  15°  is  not  common  except  near  the  equator;  more  importantly,  the  actual 
shading  error  is  so  large  for  the  combination  of  small  ft  and  high  water  turbidity  that 
measurements  should  not  be  collected  at  these  times. 


0.01  0.10  1.00 

a  (m1) 


Fig.  3.  Predictions  of  Hyper-TSRB  self-shading  in  deep  water  versus  water  absorption:  Monte 
Carlo  calculations  (triangles);  analytical  model  for  the  Hyper-TSRB  (solid  line);  and  analytical 
model  for  the  sensor  head  only  (dashed  lines).  Error  values  are  for  direct  sunlight  from  the 
solar  zenith  angles  indicated.  Analytical  values  assume  negligible  scattering;  numerical  results 
are  for  b  =  a. 

Primarily  for  academic  purposes  we  note  that  the  disagreement  between  theory  and 
numerical  result  at  small  solar  zenith  angle  in  Figs.  3  and  4  is  due  to  the  presence  of  scattering 
in  the  numerical  simulations.  The  numerical  results  converge  toward  the  theoretical  values  if 
the  amount  of  scattering  is  decreased  toward  zero.  As  demonstrated  by  Gordon  and  Ding  [1], 
scattering  does  not  have  a  large  effect  on  the  self-shading  of  a  unbuoyed  radiometer. 
However,  the  contribution  of  the  buoy  to  shading  (which  is  primarily  at  small  solar  zenith 
angles)  does  depend  significantly  on  the  amount  of  scattering  present  since  the  buoy  is 
separated  by  a  fixed  geometric  distance  away  from  the  sensor.  The  dependence  on  scattering 
of  Hyper-TSRB  self-shading  in  optically  deep  water  was  quantified  numerically  by  Leathers 
et  al.  [2].  It  is  possible  to  use  such  numerical  results  to  make  a  semianalytical  model  for  a 
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particular  instrument  that  includes  the  effects  of  scattering.  For  example,  we  can  develop  a 
semi-analytical  model  that  accounts  for  scattering  in  deep  water  as  follows.  From  Eq.  (10), 

£  =  l-exp[-(l  +  l/cos6»0w)a(z0-z,)],  (17) 

where  z0  is  given  by  Eq.  (1).  Because  the  depth  of  the  shadow  is  reduced  by  scattering,  we  can 
replace  z0  in  Eq.  (17)  with  the  reduced  shade  depth 

z'0  =rdexp(— *,  bz0),  (18) 

where  rd  is  the  radius  of  the  shading  disk  (sensor  head  or  buoy).  Shown  in  Fig.  5  is  the  self¬ 
shading  error  versus  scattering  coefficient  as  computed  with  Eqs.  (17)  and  (18)  with  k\  =  0.1. 
This  provides  a  good  model  for  the  shading  of  the  Hyper-TSRB  in  optically  deep  water; 
however,  to  apply  this  model  to  another  buoyed  instrument  it  would  be  necessary  to  compute 
a  new  value  of  k\. 


Fig.  4.  Predictions  of  Hyper-TSRB  self-shading  in  deep  water  versus  sun  position:  Monte  Carlo 
calculations  (triangles);  analytical  model  for  the  Hyper-TSRB  (solid  line);  and  analytical  model 
for  the  sensor  head  only  (dashed  lines).  Analytical  values  assume  negligible  scattering; 
numerical  results  are  for  b  -  a. 


Fig.  5.  Hyper-TSRB  self-shading  error  versus  scattering  coefficient  for  a  =  0.02  m'1.  Semi- 
analytical  values  (solid  lines)  were  obtained  with  Eqs.  (17)  and  (18)  with  kx  =  0.1,  and 
numerical  results  (asterisks)  were  computed  with  Backward  Monte  Carlo. 
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Figure  6  shows  the  Hyper-TSRB  self-shading  error  plotted  versus  total  water  depth  (zb*) 
as  computed  by  the  full  analytical  model  of  Eqs.  (13)  and  (14).  Values  are  for  Rb  =  0.2  and  a  = 
0.2  m'1  and  for  several  different  sun  positions.  Note  that  the  Hyper-TSRB  upwelling  radiance 
sensor  is  0.66  m  below  the  surface,  so  the  distance  from  the  sensor  to  the  seafloor  is  (zbot  - 
0.66).  Also  shown  in  Fig.  6  are  Hyper-TSRB  self-shading  error  results  from  Monte  Carlo 
simulations  for  b/a  =  2.  In  very  shallow  water  the  radiance  measurement  is  dominated  by  the 
component  of  light  being  reflected  off  the  seafloor.  If  the  instrument  is  very  close  to  the 
seafloor,  it  mostly  sees  its  own  shadow.  As  the  water  depth  is  increased  (i.e.,  the  instrument  is 
moved  further  from  the  bottom),  the  shading  error  in  the  radiance  reflected  off  the  bottom 
decreases,  as  does  the  overall  shading  error.  However,  as  the  water  depth  is  further  increased, 
the  water-column  component  of  the  upwelling  radiance  begins  to  dominate  and  the  overall 
shading  error  increases  toward  the  optically  deep  values.  The  existence  of  a  shading  minimum 
located  a  short  distance  from  the  seafloor  has  also  been  found  to  occur  for  upwelling 
irradiance  sensors  [15]. 


depth  (m) 

Fig.  6.  Hyper-TSRB  shading  error  as  predicted  with  the  analytical  model  compared  with  MC 
results  (triangles).  Results  shown  are  for  a  =  0.2  m1,  b  =  0.4  m'1,  Rh  =  0.2,  FOV  half-angle  = 

20°,  and  the  indicated  values  of  the  solar  zenith  angle.  The  indicated  depth  is  the  total  water 
depth;  the  sensor  is  0.66  m  below  the  sea  surface. 

The  discrepancies  in  Fig.  6  between  our  analytical  and  numerical  results  at  very  small 
water  depths  is  due  primarily  to  the  effects  of  internal  reflection  at  the  sea  surface.  When  the 
water  is  very  shallow,  some  of  the  sunlight  experiences  multiple  reflections  between  the 
seafloor  and  the  bottom  of  the  sea  surface,  making  the  overall  illumination  on  the  instrument 
less  collimated.  This  is  taken  fully  into  account  in  our  Monte  Carlo  calculations,  but  not  at  all 
in  the  analytical  results  shown  in  Fig.  6.  One  can  accommodate  internal  reflection  in  the 
analytical  model  in  the  same  way  that  atmospheric  scattering  is  accounted  for  [i.e.,  Eq.  (15)]; 
however,  this  requires  an  estimate  of  the  magnitude  of  the  deviation  of  the  in-water  light  from 
collimated. 

4.  Conclusions 

Most  commercially  available  instruments  for  measuring  in-water  upwelling  radiance  are  large 
enough  (i.e.,  with  radii  on  the  order  of  a  few  centimeters)  that  they  suffer  from  significant 
self-shading  error  over  at  least  some  portion  of  the  measured  spectrum.  The  amount  of  error 
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depends  on  instrument  dimensions,  water  optical  properties,  sun  position,  and  atmospheric 
conditions,  and  in  optically  shallow  water  also  depends  on  sensor  FOV,  water  depth,  and 
seafloor  optical  properties.  Because  radiometer  self-shading  is  minimal  when  the  water 
absorption  coefficient  is  small,  there  may  be  a  tendency  to  ignore  its  effect  in  optically  clear 
waters.  However,  even  in  the  clearest  of  waters  the  magnitude  of  the  absorption  coefficient  is 
greater  than  0.2  m'1  for  wavelengths  X  >  600  nm  and  greater  than  0.65  m'1  for  X  >  700  nm 
[12],  which  gives  deep-water  Hyper-TSRB  shading  errors  of  e  >  20%  and  e  >  50%, 
respectively,  for  6b  =  10°  and  e>  5%  and  s >  15%,  respectively,  for  #>=  30°  [Eq.  (10)]. 

The  self-shading  error  for  unbuoyed  radiometers  in  optically  deep  waters  can  be  well 
estimated  with  the  semianalytical  model  of  Gordon  and  Ding  [1].  We  have  provided  analytical 
justification  for  this  model  and  have  extended  it  to  include  the  effects  of  a  buoy  on 
instruments  such  as  the  Hyper-TSRB.  We  have  also  provided  a  semianalytical  model  for  the 
Hyper-TSRB  in  deep  water  that  takes  scattering  into  account.  For  most  practical  purposes, 
however,  the  effects  of  scattering  can  be  ignored  in  optically  deep  water. 

It  is  important  to  note  that  the  Gordon  and  Ding  model  cannot  be  applied  in  situations 
where  the  seafloor  has  a  significant  effect  on  the  measured  radiance.  We  have  developed  an 
analytical  model  for  estimating  the  self-shading  of  a  buoyed  on  unbuoyed  radiometer  close  to 
the  seafloor  [Eqs.  (11)  and  (12)]  and  linked  it  with  the  deep-water  model  with  Eq.  (13).  It 
should  be  noted  that  the  shallow-water  and  deep-water  models  do  not  bound  the  shading  error; 
Eq.  (13)  may  exhibit  an  minimum  versus  depth. 

The  error  values  provided  Eqs.  (10)  and  Eq.  (11)  are  for  a  collimated  light  incident  at  a 
particular  direction.  To  compute  the  overall  error  for  particular  illumination  conditions,  it  is 
necessary  to  properly  weigh  the  error  values  from  all  directions  [e.g.,  Eq.  (16)],  and  in  very 
shallow  water  it  may  be  necessary  to  take  internal  reflection  at  the  sea  surface  into  account. 

In  general,  the  self-shading  error  can  be  more  accurately  determined  with  Monte  Carlo 
calculations  than  with  an  analytical  model.  However,  the  analytical  approach  is  more  general 
and  is  easier  to  disseminate  and  implement.  The  disadvantage  of  the  numerical  approach  is 
that  it  requires  a  specialized  computer  code  and  that  a  large  number  of  computations  must  be 
performed  for  each  instrument.  It  is  impractical  to  expect  experimental  oceanographers  to 
write  Monte  Carlo  codes  for  this  purpose;  however,  it  may  the  reasonable  to  expect  instrument 
manufacturers  to  begin  providing  self-shading  correction  look-up  tables  or  software  with  their 
products. 

A  great  utility  of  the  analytical  model  is  its  ability  to  predict  the  amount  of  shading  that 
will  be  encountered  while  an  experiment  is  in  its  planning  stage.  Experiments  should  be 
planned  so  as  to  best  avoid  the  collection  of  upwelling  radiance  measurements  at  times  when 
self-shading  is  greatest.  For  a  given  suite  of  instrumentation,  only  the  time  of  deplojrcnent 
(e.g.,  the  position  of  the  sun)  can  be  adjusted,  but  this  may  be  sufficient.  If  an  experiment 
cannot  accommodate  the  time  restrictions,  then  it  may  be  necessary  to  invest  in  the 
development  of  smaller  sensors  or  deployment  packages  for  these  experiments.  The  analytical 
model  developed  in  this  paper  also  provides  a  tool  for  instrument  design;  the  effects  of  sensor- 
head  radius,  buoy  radius,  and  buoy  position  on  self-shading  are  all  easily  computed.  Finally  it 
should  be  noted  that  self-shading  prediction  and  correction  (by  either  analytical  model  or 
empirical  data)  requires  considerable  knowledge  of  the  environment  (e.g.,  water  properties, 
water  depth,  bottom  albedo,  and  atmospheric  conditions).  Our  analytical  model  can  be  used  to 
compute  the  sensitivity  of  self-shading  to  uncertainties  in  the  environmental  parameters  given 
a  particular  instrument’s  dimensions. 
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Abstract:  Subsurface  remote  sensing  signals,  represented  by  the  irra¬ 

diance  reflectance  and  the  remote  sensing  reflectance,  were  investigated. 
The  present  study  is  based  on  simulations  with  the  radiative  transfer 
program  Hydrolight  using  optical  properties  of  Lake  Constance  (German: 
Bodensee)  based  on  in-situ  measurements  of  the  water  constituents  and  the 
bottom  characteristics.  Analytical  equations  are  derived  for  the  irradiance 
reflectance  and  remote  sensing  reflectance  for  deep  and  shallow  water  appli¬ 
cations.  The  input  of  the  parameterization  are  the  inherent  optical  properties 
of  the  water  -  absorption  a(X)  and  backscattering  bb(X).  Additionally,  the 
solar  zenith  angle  6Sy  the  viewing  angle  By  ,  and  the  surface  wind  speed  u 
are  considered.  For  shallow  water  applications  the  bottom  albedo  Rp  and 
the  bottom  depth  zb  are  included  into  the  parameterizations.  The  result 
is  a  complete  set  of  analytical  equations  for  the  remote  sensing  signals  R 
and  Rrs  in  deep  and  shallow  waters  with  an  accuracy  better  than  4%.  In 
addition,  parameterizations  of  apparent  optical  properties  were  derived  for 
the  upward  and  downward  diffuse  attenuation  coefficients  Ku  and 
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1.  Introduction 


The  concentration  of  water  constituents  can  be  derived  by  optical  remote  sensing  making  use 
of  the  spectral  shape  of  the  reflected  sunlight.  The  models  of  [1],  [2],  or  [3]  for  the  irradiance 
reflectance  and  remote  sensing  reflectance,  R  and  Rrs ,  use  as  parameters  the  absorption  a(X) 
and  backscattering  coefficient  bh(X): 


Roo 


r 


H 

a  +  bb 


=  fx 


f 


■1  bb 
a  +  bh 


(1) 

(2) 
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The  equations  are  valid  just  below  the  water  surface.  The  higher  the  backscattering,  the  higher 
the  reflectance  and  the  higher  the  absorption,  the  lower  the  reflectance.  f°  is  the  proportionality 
factor  for  the  irradiance  reflectance  and  /f  for  the  remote  sensing  reflectance  of  an  infinitely 
deep  water  body  (index  «>).  The  g-factor  is  defined  as  the  ratio  of  the  upwelling  irradiance 
and  the  upwelling  radiance:  Q  =  j*.  For  an  isotropic  upwelling  radiance  distribution,  Q  =  n  sr. 
Thus,  the  g-factor  is  a  measure  of  the  anisotropy  of  the  light  field  distribution.  In  natural  waters, 
g  is  typically  around  5  sr.  Using  the  definitions  of  the  irradiance  reflectance  and  remote  sensing 
reflectance,  R  =  and  Rrs  =  jps  the  g-factor  is  also  the  ratio  of  the  irradiance  reflectance  to 
the  remote  sensing  reflectance  according  to  Eq.  (2). 

These  approximations  are  valid  for  infinitely  deep  water,  where  only  the  water  body  con¬ 
tributes  to  the  reflected  signal.  For  open  ocean  case-1  waters,  constant  proportionality  factors 
are  commonly  used  and  are  sufficient  for  many  applications,  for  example  f°  =  0.33  [1]  and 
/T  —  0.095  [4].  Due  to  the  influences  of  the  sun  and  surface  geometry  on  the  reflectances,  a 
parameterization  of  the  factors  can  be  developed  including  these  aspects  and  the  inherent  op¬ 
tical  properties  as  well  [5,  6].  But  in  case-2  waters,  the  factors  f°  and  are  not  constant  and 
can  vary  in  time  and  location  [7].  [8]  used  for  his  determination  at  Lake  Constance  (German: 
Bodensee)  a  combination  of  the  model  of  [5]  and  [9]. 

In  addition,  if  remote  sensing  data  are  analyzed  including  optically  shallow  waters,  the  bot¬ 
tom  depth  zb  and  the  bottom  albedo  Rb  have  to  be  taken  into  account  [10].  Following  [11], 
different  authors  including  [12],  [13],  or  [14]  have  formulated  approximations  of  the  irradiance 
reflectance  for  shallow  water.  Their  equations  result  from  a  two-flow  irradiance  approximation 
including  the  bottom  influences: 

R  =  Roo  (l  -  e-2*z*)  +  Rb£~2Kzr 

K  is  the  diffuse  attenuation  coefficient  of  the  water  column  and  is  equal  for  the  downward  and 
upward  directions.  According  to  [15],  this  equation  can  be  transformed  to  the  remote  sensing 
reflectance  Rrs  as  follows: 

However,  in  reality  the  diffuse  attenuation  coefficients  of  the  upwelling  and  downwelling  light 
are  not  equal.  To  get  a  more  accurate  expression  the  effective  attenuation  coefficient  is  divided 
into  an  upwelling  and  a  downwelling  part.  The  upwelling  part  distinguishes  between  radiation 
reflected  in  the  water  column  (index  W)  and  from  the  bottom  (index  B).  This  results  in  the 
following  equations,  as  mentioned  in  [15]  and  [16]: 

R  =  R„  (l  -  e~(*''+*“  tv)zs)  +  Rge'-^+x^B 

R„  =  *».-[' -««p{-(*j  +  ^|)a}] 

The  viewing  angle  just  below  the  water  surface,  9Vy  is  included  in  Eq.  (4)  in  the  upward  attenu¬ 
ation  due  to  the  dependence  of  the  upwelling  radiance  on  the  viewing  position. 

Based  on  these  equations  new  analytical  parameterizations  were  developed  for  the  re¬ 
flectances  of  deep  and  shallow  water,  the  upwelling,  and  the  downwelling  attenuation  coef¬ 
ficients  to  obtain  an  invertable  equation  for  remote  sensing  data.  Simulations  were  made  with 
the  well-established  and  validated  radiative  transfer  program  Hydrolight  (version  3.1)  [17,  18] 
for  case-2  waters  and  the  results  fitted  to  Eq.  (3)  and  (4).  The  program  code  was  optimized  for 


(3) 

(4) 
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the  test  site  Lake  Constance,  located  in  southern  Germany  at  the  borders  with  Switzerland  and 
Austria.  The  optical  properties  of  Lake  Constance  were  investigated  by  [19]  and  [8]  and  were 
included  into  the  source  code  of  Hydrolight.  The  simulations  with  Hydrolight  were  performed 
not  only  over  the  natural  range  of  the  concentrations  of  the  water  constituents  found  for  Lake 
Constance,  but  also  below  and  above  this  range  to  cover  a  more  general  range  of  concentra¬ 
tions.  This  extends  the  validity  of  the  developed  parameter! zations  for  a  wide  number  of  case-2 
waters. 

2.  Radiative  transfer  model 

For  the  forward  simulation  of  the  underwater  light  field  the  radiative  transfer  program  Hydro¬ 
light  (version  3.1)  washed.  The  model  is  explained  in  detail  by  [17]  and  [18].  The  optical 
properties  are  selected  for  case-2  waters,  such  as  the  test  site  at  Lake  Constance,  and  are  given 
for  three  kinds  of  water  constituents:  phytoplankton  (index  P),  suspended  matter  (index  X),  and 
gelbstoff  (index  Y).  Gelbstoff  is  also  known  as  Colored  Dissolved  Organic  Matter  (CDOM), 
yellow  substance,  or  gilvin.  Thus,  the  total  absorption  and  backscattering  coefficients,  a(A)  and 
bf,( A),  are  the  sums  of  all  contributions  of  each  water  constituted  and  pure  water  (index  W) 
itself: 

tf(A)  =  «w(A)+tfp(A)  +  ax(A)+0K(A) 
bb  ( A )  =  bby  w  ( A )  +  A )  +  bbyx  (A )  +  bby  (A ) 

Due  to  the  symmetric  scattering  properties  of  water,  the  backscattering  coefficient  of  pure  wa¬ 
ter  can  be  determined  from  the  scattering  coefficient  bw :  bby/  =  jbw-  The  absorption  aw( A) 
and  scattering  coefficient  bw( A)  of  pure  water  are  taken  from  [20].  The  absorption  and  back- 
scattering  coefficients  of  phytoplankton  and  suspended  matter  are  the  product  of  the  specific 
absorption  and  backscattering  coefficients  and  the  concentrations.  The  specific  absorption  of 
phytoplankton  aP( A)  is  given  after  [8]  by  the  mean  value  of  measurements  of  [19]  at  Lake 
Constance.  Absorption  of  suspended  matter  was  investigated  for  Lake  Constance  by  [8].  He 
found  no  specific  absorption  and  therefore  ax  is  set  to  zero.  The  absorption  of  gelbstoff  is  ap¬ 
proximated  by  an  exponential  function  [21]  with  an  exponent  S  =  0.014  nm~*  at  a  reference 
wavelength  Ao  of  440  nm  after  [22].  The  scattering  in  the  water  is  mainly  driven  by  the  amount 
of  suspended  matter.  The  influence  of  the  particulate  fraction  of  phytoplankton  on  the  scattering 
is  included  in  the  value  of  the  scattering  and  backscattering  coefficients  of  suspended  matter  as 
investigated  by  [8]  and  therefore,  bbyp  =  0.  There  was  also  no  dependence  of  the  scattering  and 
backscattering  coefficients  found  at  the  wavelengths  between  400  and  800  nm.  Hence,  the  con¬ 
stant  value  of  the  specific  backscattering  coefficient  of  suspended  matter  b*b  X  =  0.0086  m2/g 
obtained  for  Lake  Constance  was  used.  The  ratio  of  the  specific  scattering  to  backscattering 
coefficient  is  0.019  [8],  which  is  the  same  as  found  by  [23]  in  San  Diego  harbor.  Therefore,  this 
phase  function  was  chosen  for  all  simulations.  Gelbstoff  is  assumed  not  to  scatter  light  because 
its  pigments  are  totally  dissolved  in  the  water.  Finally,  the  following  equations  are  obtained  for 
the  absorption  and  the  backscattering  coefficients: 

a(  A)  =  flvv(A)  +aP(X)Cp  4-  ar(Ao)e~5'*-^  (5) 

b„(X)  =  l-bw(X)  +  b'hJ(Cx  (6) 

Cp  is  the  concentration  of  phytoplankton,  which  is  given  as  the  sum  of  chlorophyll-a  and 
phaeophytin-a  concentration  in  units  of  fig/ 1.  Cx  is  the  concentration  of  suspended  matter  in 
units  of  mg/1. 

The  impact  of  inelastic  processes  in  natural  water  was  recently  investigated  by  [24]  and  found 
to  contribute  significantly  to  the  irradiance  reflectance  and  remote  sensing  reflectance.  Hence, 
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the  fluorescence  of  chlorophyll  and  gelbstoff  as  well  as  Raman  scattering  were  included  in  all 
simulations  using  the  default  efficiencies  of  Hydrolight  (version  3.1).  The  quantum  efficiency 
of  chlorophyll  fluorescence  was  set  to  2%  and  the  emission  function  was  approximated  by  a 
Gaussian  function  centered  at  the  wavelength  685  nm  with  10.6  nm  full  width  at  half  maxi¬ 
mum.  The  fluorescence  of  gelbstoff  was  described  by  the  spectral  fluorescence  quantum  effi¬ 
ciency  function  defined  by  [25]  between  310  and  490  nm.  The  quantum  efficiency  took  values 
of  about  0.9  to  1 .9%  in  this  wavelength  interval.  The  Raman  scattering  was  approximated  using 
a  wavenumber  redistribution  function  expressed  by  the  sum  of  four  Gaussian  functions  repre¬ 
senting  the  number  of  shifts  for  a  scattered  photon.  For  more  details  see  [18]. 

The  water  constituents  are  in  general  not  homogeneously  distributed  with  depth  in  natural 
waters.  Thus,  to  match  on  average  the  real  situations  of  the  test  site  Lake  Constance,  depth 
profiles  of  the  water  constituents  were  included  in  all  the  simulations.  More  than  500  depth 
profiles  were  measured  at  Lake  Constance  between  1986  and  1996.  The  data  were  analyzed 
and  mean  profiles  were  determined  for  phytoplankton  and  suspended  matter  [8,  26].  The  mini¬ 
mum  concentration  of  phytoplankton  was  1.0  /ig/1  used  for  determining  the  depth  profile.  The 
dependence  of  the  concentration  on  the  depth  z  can  be  expressed  as 

C(z)  ^Q+C^exp  j-|  (7) 

where  C{z)  is  the  concentration  of  phytoplankton  or  suspended  matter.  For  gelbstoff,  no  depth 
dependence  was  found.  The  values  of  the  coefficients  Co,  zmax .  cr>  and  n  are  listed  in  table 
1.  If,  for  example,  a  constant  value  of  the  concentration  for  all  depths  is  used,  the  irradiance 
reflectance  is  underestimated  by  12  to  15%  for  concentrations  of  2  to  5  /xg/1  phytoplankton  and 
2  to  5  mg/I  suspended  matter. 


Table  1 .  Values  of  the  constant  factors  of  phytoplankton  and  suspended  matter  for  the  depth 
profile  in  Eq.  (7). 


Co 

Zmax  (m) 

<r(m) 

n 

Phytoplankton 

1.0  /ig/1 

2.9 

9.6 

3.0 

Suspended  matter 

0.9  mg/1 

12.9 

10.7 

2.3 

For  the  simulations  from  400  to  750  nm  with  steps  of  5  nm,  the  concentrations  of  the  water 
constituents  at  the  surface  were  varied  beyond  their  natural  ranges  at  Lake  Constance  (see 
figures  1  and  2).  The  values  are  given  in  table  2.  To  get  a  suitable  depth  profile  for  the  lowest 
concentrations  of  phytoplankton  and  suspended  matter  the  coefficient  Co  in  table  1  was  adjusted 
to  0.5  /ig/1  and  0.5  mg/1  respectively. 


Table  2.  Concentrations  of  the  water  constituents  for  the  simulations  with  Hydrolight. 


Cp  (/ig/1) 

0.5 

20.0 

1.0 

40.0 

1.5 

60.0 

2.0 

80.0 

3.0 

100.0 

5.0 

7.0 

10.0 

Cx  (mg/I) 

0.5 

1.0 

1.5 

2.0 

3.0 

4.0 

5.0 

6.0 

7.0 

9.0 

10.0 

30.0 

50.0 

artAoXm-1) 

0.05 

0.10 

0.15 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 

0.80 

0.90 

1.00 

1.30 

2.50 

4.00 

5.00 
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Fig.  1.  Distribution  of  the  concentration  of  suspended  matter  against  phytoplankton  for  the 
simulations  with  Hydrolight. 


For  the  bottom  albedo,  spectra  of  sand  and  green  algae  included  in  Hydrolight  were  used. 
Additionally,  bottom  albedo  spectra  were  measured  at  Lake  Constance  for  sediment  and  macro¬ 
phyte  covered  beds  with  the  Hydrological  Spectral  Radiometer  HYDRA  [27].  The  sediment 
was  gray  colored  mud  consisting  of  about  50%  inorganic  particles  and  50%  organic  matter  typ¬ 
ically  for  Lake  Constance.  The  average  particle  size  is  around  10  /im.  The  reflectance  of  the 
macrophyte  was  measured  above  a  patch  of  the  species  Potamogeton  Pectinatus  L.  with  senes¬ 
cent  leaves.  The  bottom  reflectances  used  for  the  simulations  are  shown  in  figure  3.  Bottom 
depth  was  set  to  1,3,  6,  10,  20,  30  m,  and  infinitely  deep  water.  For  simulating  the  sun  and 
sky  conditions  the  model  of  [28]  included  in  Hydrolight  was  used.  Clear  sky  conditions  were 
chosen  with  varying  subsurface  sun  zenith  angle  ft.  The  angles  were  8,  14,  21,  27,  34,  39,  43, 
and  46°.  The  forward  simulations  were  performed  for  values  of  the  surface  wind  speed  u  of  0, 
5,  and  10  m/s  using  the  sea  surface  statistical  model  incorporated  into  Hydrolight.  This  surface 
representation  is  based  on  the  wave-slope  wind-speed  laws  of  [29,  30]  and  thus  includes  both 
gravity  and  capillary  wave  effects. 

Altogether,  over  1400  spectra  were  simulated.  The  output  of  the  Hydrolight  simulations  is 
given  for  different  subsurface  viewing  angles  ft  ranging  between  8°  and  46°  in  the  same  manner 
as  the  sun  zenith  angle. 

3.  Results 

The  output  of  Hydrolight  of  the  irradiance  reflectance  and  remote  sensing  reflectance,  R  and 
Rrs ,  were  investigated  to  find  parameterizations  for  the  inherent  optical  properties  as  well  as 
for  the  sun  and  viewing  geometry  using  Eq.  (3)  and  (4).  The  unknown  variables  in  these  equa¬ 
tions  are  the  irradiance  reflectance  and  remote  sensing  reflectance  of  the  water  body,  /?«*>  and 
Rrs% oo,  and  the  diffuse  attenuation  coefficients  Kj,  Ku%w,  and  KUjb .  For  each  of  these  unknown 
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Fig.  2.  Distribution  of  the  concentration  values  of  gelbstoft  (top),  suspended  matter  (center), 
and  phytoplankton  (bottom)  for  the  simulations  with  Hydrolight. 


factors  new  parameterizations  were  developed  based  on  the  inherent  optical  properties,  the  so¬ 
lar  and  viewing  angles,  and  the  surface  wind  speed.  Wavelength  dependence  was  included  in 
the  inherent  optical  properties.  For  all  simulations  chlorophyll  and  gelbstoff  fluorescence  was 
considered  as  mentioned  above.  For  the  determination  of  the  parameterizations,  wavelengths 
around  the  strongly  peaked  fluorescence  emission  of  chlorophyll  from  660  to  715  nm  were  ex¬ 
cluded  to  give  a  better  spectral  fit.  The  parameterizations  were  generated  by  fitting  the  simulated 
values  using  the  Levenberg-Marquardt  multivariate  optimization  technique  yielding  regression 
coefficients.  These  regression  coefficients  allow  calculation  of  the^ reflectances  and  attenuation 
coefficients  using  analytical  equations.  The  mean  relative  error  8  documents  the  accuracy  of 
the  analytical  parameterizations  and  was  estimated  by 

-_X-Xo  AX 
X0  X0  ’ 

where  X  represents  the  value  calculated  by  the  parameterization  and  Xq  the  value  of  the  simu¬ 
lation. 
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Fig.  3.  Bottom  reflectance  spectra  used  for  the  forward  simulations  in  Hydrolight. 


3.1 .  Deepwater 

3.1.1.  Irradiance  reflectance 

The  underwater  irradiance  reflectance  for  deep  water  /?«,  at  depth  z  =  0  is  parameterized  using 
Eq.  (1).  The  factor  f°  was  analyzed  for  its  dependence  on  x ,  w,  and  0S.  The  simulations  yield  a 
non-linear  dependence  on  the  factor  x  (Fig.  4)  and  the  subsurface  solar  zenith  angle  0S  (Fig.  5 
right)  and  a  linear  dependence  on  the  surface  wind  u  (Fig.  5  left).  For  the  investigation  of  the 
dependence  on  surface  wind,  additional  calculations  were  made  for  wind  speed  values  ranging 
from  0  to  30  m/s  in  steps  of  1  m/s.  The  following  parameterization  was  found  to  be  suitable: 


=  f(x,et,u)x  =  f(x)f(6s)f(u)x 


(8) 


Using  fewer  coefficients  for  the  factor  x  results  in  a  significantly  lower  correlation.  The  coeffi¬ 
cients  p\  to  p(y  were  determined  using  N  =  22184  model  results.  The  results  of  the  regression 
are  listed  in  table  3.  The  errors  of  using  a  constant  factor  f°  is  illustrated  by  figure  4.  The  dashed 
line  corresponds  to  the  value  of  f°  =  0.33  by  [  1].  For  high  values  of  xy  which  means  high  back- 
scattering  or  high  concentration  of  suspended  matter,  the  error  increases  up  to  100%.  With  the 
new  parameterization  the  error  is  reduced  significantly.  Figure  6  shows,  on  the  left  hand  side, 
the  calculated  plotted  against  the  simulated  values.  The  black  crosses  are  the  estimated  values 
with  Eq.  (8)  and  the  blue  points  are  the  values  derived  by  the  previous  model  of  [8].  The  dis¬ 
tribution  of  the  relative  error  8  between  the  simulated  and  predicted  values  of  with  Eq.  (8) 
(Fig.  6  right)  shows  a  normal  distribution  with  a  mean  value  of  0.04  while  the  mean  error  using 
a  constant  f°  =  0.33  is  -0.25.  The  new  parameterization  gives  also  much  better  results  than 
models  of  [5]  and  [6]  which  include  the  sun  zenith  angle  to  estimate  the  factor  f°  (graph  not 
shown).  The  previous  model  of  [8]  for  Lake  Constance  results  in  a  mean  relative  error  of  0.08 
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(Fig.  6  right).  But  for  irradiance  reflectances  greater  than  20%  the  mean  relative  error  is  0.22. 
For  these  situations  of  high  backscattering  due  to  high  amount  of  suspended  matter,  the  new 
parameterization  results  in  a  mean  relative  error  below  1%. 


Fig.  4.  Irradiance  reflectance  for  infinitely  deep  water  simulated  with  Hydrolight  (N  = 
22184)  depending  on  x  =  5^  and  the  approximation  of  [1 )  for  a  factor  f°  =  0.33  (dashed 
line). 

The  advantage  of  the  new  parameterization  is  the  separation  of  the  dependences  on  the  in¬ 
herent  optical  properties  and  the  sun  and  surface  geometry.  This  allows  the  influence  of  the 
variables  on  the  remote  sensing  signal  to  be  analyzed  separately.  The  surface  wind  speed  has 
the  weakest  influence.  If  it  is  neglected,  the  error  is  below  1%.  The  influence  of  the  sun  position 
is  greater:  the  variation  of  the  irradiance  reflectance  is  about  15%  for  a  subsurface  solar  zenith 
angle  from  0  to  25°  and  about  30%  for  0  to  40°. 

3.1.2.  Remote  sensing  reflectance 

The  determination  of  the  remote  sensing  reflectance  from  the  irradiance  reflectance  is  possible 
using  the  Q-factor,  Rrs  oo  =  The  (Mactor  is  determined  by  the  angular  distribution  of  the 
light  field  under  water.  Therefore,  a  parameterization  of  Q  =  Q(6S}  0yyu)  seemed  to  be  suitable. 
All  data  points  were  analyzed,  but  no  suitable  parameterization  was  found.  The  reason  is  that  the 
angular  distribution  of  Q  is  controlled  also  significantly  by  the  inherent  optical  properties  and 
their  concentrations:  Q  =  Q(0s,0v,uyx).  Thus,  an  equation  for  the  remote  sensing  reflectance 
in  deep  water  was  established  that  is  similar  to  the  equation  for  the  irradiance  reflectance,  but 
with  different  values  of  the  coefficients.  The  factor  /T  is  derived  as  a  function  of  separated 
variables.  In  addition  to  the  dependences  on  x,  0S>  and  w  the  remote  sensing  reflectance  varies 
with  the  subsurface  viewing  angle  6V.  Simulations  using  different  values  of  0V  are  shown  in 
Fig.  7.  The  data  points  can  be  fitted  with  a  function  proportional  to  The  variation  of  the 
remote  sensing  reflectance  is  about  10%  for  a  subsurface  viewing  angle  from  0  to  25°.  This  is 
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Fig.  5.  Dependence  of  the  irradiance  reflectance  for  infinitely  deep  water  on  surface  wind 
(left)  and  subsurface  solar  zenith  angle  (right).  The  concentrations  of  the  water  constituents 
are  C/>  =  3  Jig/1,  C\  =  3  mg/1,  and  fly(Ao)  =  0.2m_1 . 
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Fig.  6.  Left:  irradiance  reflectance  calculated  by  Eq.  (8)  (black  crosses)  and  by  the  model 
of  [8]  (blue  points)  against  the  simulated  values  for  infinitely  deep  water.  The  1:1  line  is 
plotted  in  red.  Right:  distribution  of  the  relative  errors  for  the  approximation  of  [1]  (white 
bars),  (8)  (gray  bars),  and  of  the  new  parameterization  of  Eq.  (8)  (cross  hatched  bars). 


accounted  for  using  the  following  parameterization: 

=  fHx,ds,uA)x  =  fHx)fH0s)f\u)fHev)x 

=  Pl  (1  +  P2X  +  P1X2  +  P4X3) 

*(1+,W)*  <9> 

The  results  of  the  regression  are  listed  in  table  3.  N  =  177472  model  results  were  used  to 
calculate  the  coefficients  of  the  equation.  Figure  8  shows  the  calculated  values  plotted  against 
the  simulated.  The  mean  relative  error  is  about  0.02. 
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Fig.  7.  Dependence  of  the  remote  sensing  reflectance  for  infinitely  deep  water  on  the  sub¬ 
surface  viewing  angle  ft,.  The  concentrations  of  the  water  constituents  are  Cp  =  3  jzg/1, 
Cx  -  3  mg/1,  and  ay(Ao)  —  0.2m”1. 


Table  3.  Coefficients  for  the  irradiance  reflectance  of  deep  water  for  Eq.  (8)  and  for  the 
remote  sensing  reflectance  of  deep  water  for  Eq.  (9). 


Ro.  of  Eq.  (8) 

/?rj,»(sr l)  of  Eq.  (9) 

P\ 

0.1034  ±0.0014 

0.0512  ±  0.0001  sr”1 

P2 

3.3586  ±  0.0305 

4.6659  ±0.0174 

P3 

-6.5358  ±  0.0808 

-7.8387  ±  0.0434 

PA 

4.6638  ±  0.0649 

5.4571  ±0.0345 

PS 

2.4121  ±0.0443 

0.1098  ±0.0018 

P6  (s/m) 

-0.0005  ±  0.0001 

-0.0044  ±  0.0000 

Pi 

- 

0.4021  ±  0.0020 

3.2.  Diffuse  attenuation  coefficient 

The  reflectances  of  deep  water  are  the  input  for  the  shallow  water  Eq.  (3)  and  (4).  To  employ 
these  equations  it  is  necessary  to  estimate  the  diffuse  attenuation  coefficients.  The  diffuse  at¬ 
tenuation  describes  the  loss  of  up-  and  downwelling  irradiance  within  a  thin  layer  in  the  water. 
This  loss  depends  on  absorption,  scattering,  and  the  isotropy  of  the  light  field.  The  latter  is 
parameterized  by  the  mean  cosine  Ji.  For  a  totally  isotropic  light  distribution  Ji  is  0  and  for  a 
collimated  beam  in  direction  0  the  mean  cosine  has  the  value  JI  =  cos  6  (see  [  1 8]  for  example). 
In  clear  sky  conditions,  just  below  the  water  surface  the  distribution  of  the  light  is  mainly  af¬ 
fected  by  the  direct  beam  of  the  sun.  Thus,  the  mean  cosine  is  approximately  the  cosine  of  the 
subsurface  solar  zenith  angle  6S  in  the  upper  water  layers. 
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Fig.  8.  Left:  remote  sensing  reflectance  calculated  by  Eq.  (9)  against  the  simulated  values 
for  infinitely  deep  water.  The  1:1  line  is  plotted  in  red.  Right:  distribution  of  the  relative 
errors. 


3.2. 1 .  Downward  diffuse  attenuation 

Kj  depends  on  the  absorption  and  backscattering  as  well  as  on  the  solar  zenith  angle  as  shown 
in  figure  9  (left).  After  [31],  the  parameterization  for  Kd  is 


Kd  =  Ko 


a  +  h 

cos  6S 


(10) 


The  simulated  values  of  Kd  range  from  about  0.1  m”1  to  10.6  m~l  with  a  mean  value  of  0.7 
m“ 1 .  The  regression  oiN  =  72558  data  points  gives  a  value  of  Kq  =  1 .0546  ±  0.000 1 .  The  mean 
relative  error  is  5  =  -0.01 .  The  distribution  of  the  relative  errors  is  shown  in  figure  9  (right). 


Fig.  9.  Downward  diffuse  attenuation  coefficient  of  72558  simulations  with  Hydrolight. 
Left:  dependency  on  Right:  distribution  of  the  relative  errors  between  calculated  and 
simulated  values. 


3.2.2.  Upward  diffuse  attenuation 

The  investigation  on  the  upwelling  diffuse  attenuation  coefficient  is  first  done  for  the  infinitely 
deep  water  body  to  find  the  best  parameterization.  Figure  10  shows  the  dependence  of  Ku  on 
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absorption  and  backscattering  (left)  and  on  the  subsurface  solar  zenith  angle  (right).  The  graph 
on  the  right  hand  side  indicates  the  following  dependence  of  Ku  on  the  subsurface  solar  zenith 
angle  6S:  Ku  «=  On  the  left  hand  side  the  dependence  on  absorption  and  backscattering  is 
plotted.  The  colors  indicate  the  concentration  of  suspended  matter  Cx  which  is  directly  linked 
to  the  backscattering  coefficient  as  described  in  Eq.  (6).  Generally,  a  linear  dependence  on  the 
sum  of  absorption  and  backscattering  can  be  assumed:  Ku  «  (0  +  6,,).  However,  for  concentra¬ 
tions  of  suspended  matter  of  Cx  <  3  0  mg/1  the  relationship  differs  increasingly  from  a  linear 
dependence.  The  upward  diffuse  attenuation  coefficient  takes  higher  values  for  a  lower  amount 
of  scattering  particles  in  the  water.  This  is  because  few  photons  are  scattered  upwards,  result¬ 
ing  in  an  anisotropic  light  field  [32].  To  correct  for  this  effect  an  additional  term  is  included 
depending  on  jc  =  The  following  equation  for  the  upward  diffuse  attenuation  coefficient 
shows  the  best  fit  when  used  in  Eq.  (3)  and  (4). 


«r.-(.+6.)(i  «)'■(.  +  «^)  <ii) 

For  the  simulations  of  an  infinitely  deep  water  body  the  mean  relative  error  was  8  =  0.13  for 
N  =  22184  points.  To  separate  the  influence  of  photons  reflected  by  the  water  column  and 
the  bottom,  two  different  upward  diffuse  attenuation  coefficients,  Ku>w  and  Kuj *,  are  used  for 
shallow  water  applications.  Thus,  four  coefficients  and  with  i—  1,2  were  determined 
by  fitting  the  entire  Eq.  (3)  and  (4).  The  results  of  the  regression  are  given  in  table  4.  Since  the 
output  of  Hydrolight  is  the  total  upward  diffuse  attenuation  coefficient,  which  is  not  the  sum  of 
Ku,w  and  no  mean  relative  errors  can  be  specified. 


Fig.  10.  Dependency  of  the  upward  diffuse  attenuation  coefficient  on  the  sum  of  absorption 
and  backscattering  (left)  and  subsurface  solar  zenith  angle  (right).  The  points  on  the  left  are 
for  0S  =  8°  with  colors  representing  the  concentration  of  suspended  matter  and  the  curve 
on  the  right  is  for  C>  =  1  /xg/l,  Cx  =  1  mg/1,  and  ay  (Ao)  =  0.2m”1. 


3.3.  Shallow  water 

Putting  all  the  above  results  together  the  shallow  water  reflectances  can  be  calculated  using 
Eq.  (3)  and  (4).  Additional  coefficients  A\  and  A  2  are  introduced  to  adapt  the  equations  to  the 
simulated  situations. 
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Fig.  1 1.  lrradiance  reflectance  of  shallow  water.  Comparison  of  simulated  values  Ro  and 
estimated  values  R  (left)  with  the  1:1  line  in  red:  the  green  points  are  the  values  for  wave¬ 
lengths  from  660  to  715  nm.  Distribution  of  the  relative  error  between  simulated  and  esti¬ 
mated  irradiance  reflectances. 


3.3. 1 .  lrradiance  reflectance 

With  the  new  parameterizations  of  the  diffuse  attenuation  coefficients  and  the  factors  A\  and 
Aiy  the  irradiance  reflectance  can  be  expressed  by 

*  -  R-[i-','“|,{'(^+(l+')r"(,+S))(“+',‘)a}] 

+  *A“p{-(^;  +  (i+»l''(it^))(*+«»}  02) 

The  values  of  A\  and  A  2  were  determined  by  fitting  N  =  72558  simulated  data  points.  The 
resulting  coefficients  are  listed  in  Table  4.  In  Fig.  1 1  the  estimated  irradiance  reflectance,  /?.  is 
plotted  against  the  simulated  reflectance,  Ro,  for  all  cases.  The  distribution  of  the  relative  error 
-  with  a  mean  error  of  5  =  0.02  -  is  also  shown  in  Fig.  1 1 .  For  comparison,  the  original  Eq.  (3) 
of  [13]  using  Eq.  (8)  for  the  irradiance  reflectance  of  the  water  column  gives  a  relative  mean 
error  of  5  =  0.06. 

The  green  points  in  Figure  1 1  are  the  values  for  wavelengths  from  660  to  715  nm.  Although 
these  wavelengths  were  excluded  for  algorithm  development  due  to  the  influence  of  chlorophyll 
fluorescence,  the  estimation  with  Eq.  (12)  using  the  parameters  of  Table  4  is  in  fair  agreement 
with  the  simulated  values.  The  mean  relative  error  for  these  wavelengths  is  5  =  —0.12.  This 
means  that  the  new  parameterization  can  be  applied  also  to  these  wavelengths  with  an  underes¬ 
timation  of  about  12%. 

The  spectral  shape  of  three  examples  is  shown  in  Fig.  12  in  a  range  from  400  to  750  nm 
with  the  relative  errors  compared  to  the  simulations  of  Hydrolight.  The  numbers  in  the  figure 
correspond  to  the  following  conditions: 

1.  Bottom  type  is  sediment  at  a  depth  of  zb  =  5  m;  the  concentration  of  the  water  con¬ 
stituents  are  phytoplankton  C/>  =  10.8  /ig/1,  suspended  matter  Cx  =  50.0  mg/1,  and  gelb- 
stoflf  ay  (440nm)  =  0.2  m_1 ;  the  subsurface  solar  zenith  angle  is  6S  =  27°;  the  wind  speed 
is  u  =  1  m/s. 

2.  Macrophytes  at  zb  =  6  m;  C/>  =  2.5  jig/1,  Cx  =  7.0  mg/l,  and  ay(440nm)  =  0.3  m_1; 
0S  =  33°;  u  —  0  m/s. 
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3.  Sediment  at  zb  =  5  m;  Cp  =  1.0  jUg/1,  C*  =  1.0  mg/l,  and  0y  (440nm)  =  0.05  m”1;  ft  = 
27°;  m  =  1  m/s. 

The  agreement  between  simulation  and  calculation  with  Eq.  (12)  is  very  good  in  each  case.  The 
relative  error  (Fig.  12:  right)  is  below  5%  over  the  entire  spectral  range,  except  for  wavelengths 
around  685  nm.  This  is  due  to  the  fluorescence  of  chlorophyll  which  is  not  parameterized  in 
these  analytical  equations.  The  differences  in  the  other  parts  of  the  spectra  are  mainly  caused  by 
the  fluorescence  of  gelbstoff  which  affects  mostly  the  green  part  of  the  visible  spectrum  [24]. 


Fig.  12.  Irradiance  reflectance  of  shallow  water  for  the  spectral  range  from  400  to  750 
nm  for  three  different  cases.  Left:  comparison  of  simulated  (dotted  lines)  and  estimated 
values  (solid  lines);  right:  relative  errors.  The  numbers  refer  to  the  following  situations:  (1) 
sediment  at  zb  =  5  m,  Cp  =  10.8  /zg/1,  Cx  =  50.0  mg/l,  ay(440nm)  =  0.2  m~l,  ft  —  27°, 
u  =  1  m/s;  (2)  macrophytes  at  zb  =  6  m,  Cp  =  2.5  jig/1,  Cx  =  7.0  mg/l,  ay(440nm)  — 
0.3  m"1,  ft  =  33°,  u  =  0  m/s;  (3)  sediment  at  zb  =  5  m,  Cp  =  1.0  jig/1,  Cx  =  10  mg/l, 
ay  (440nm)  =  0.05  m"1 ,  ft  =  27°,  u  =  1  m/s. 


3.3.2.  Remote  sensing  reflectance 

The  remote  sensing  reflectance  Rrs  can  be  expressed  with  a  similar  approximation  as  the  irra- 
diance  reflectance,  but  with  an  additional  dependence  on  the  subsurface  viewing  angle  ft,. 


=  Rn 


(V2!£  +  (,  +,)«»  (,  +  JStY)  £+£„)] 

(  \  cos  ft  \  cos  ft  /  /  cos  ft  J  J 


4-  A 


cosft,  xr,  o 
«‘^+(i+x> 


< 13> 

\  cos  ft  /  /  cos  ft,  J 


The  results  of  the  regression  analysis  of/V  =  580464  numbers  of  observations  are  listed  in  Table 

4  below.  Figure  13  shows  the  estimated  against  the  simulated  values  and  the  relative  error,  with 
a  mean  error  of  S  =  0.03.  For  comparison,  the  equation  of  [15]  gives  a  relative  mean  error  of 
-0.09. 

As  mentioned  for  the  case  of  the  irradiance  reflectance,  the  green  points  in  Fig.  13  are  the 
values  for  wavelengths  from  660  to  715  nm.  The  correlation  between  the  estimated  and  sim¬ 
ulated  values  are  very  good  here  as  well.  The  mean  relative  error  for  these  wavelengths  is 

5  =  —0.13.  This  means  that  the  new  parameterization  can  be  applied  also  to  these  wavelengths 
with  an  underestimation  of  about  13%. 
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Table  4.  Coefficients  for  the  irradiance  and  remote  sensing  reflectance  of  Eq.  (12)  and  (13) 
for  shallow  water. 


KofEq.  (12) 

A\ 

1.0546  ±  0.0038 

*b 

1.0546  ±  0.0001 

K\,w 

1.9991  ±0.0305 

*2 

0.2995  ±  0.0122 

M 

0.9755  ±  0.0013 

K'l.fi 

1.2441  ±0.0209 

K2  ,B 

0.5182  ±  0.0036 

W„(sr~‘)  of  Eq.  (13) 
1.1576  ±0.0014 
1.0546  ±  0.0001 
3.5421  ±0.0152 
-0.2786  ±  0.0030 
1.0389  ±  0.0004  Sr"1 
2.2658  ±  0.0076 
0.0577  ±  0.0009 


Fig.  13.  Remote  sensing  reflectance  of  shallow  water.  Comparison  of  simulated  values  Rrs$ 
and  estimated  values  Rrs  (left)  with  the  1:1  line  in  red;  the  green  points  are  the  values  for 
wavelengths  from  660  to  715  nm.  Distribution  of  the  relative  error  between  simulated  and 
estimated  remote  sensing  reflectances. 


The  spectral  shape  of  the  remote  sensing  reflectance  is  shown  in  Fig.  14  in  the  same  way  as 
explained  for  the  irradiance  reflectance.  Two  graphs  are  included  in  the  figure,  one  for  a  subsur¬ 
face  viewing  angle  of  0V  =  7°  and  one  for  Gv  =  27°.  The  calculated  values  of  the  remote  sensing 
reflectance  agree  very  well  with  the  simulations.  The  relative  error  is  below  5%  except  for  situ¬ 
ation  number  3,  where  the  error  is  about  10%  around  600  nm.  The  reason  for  the  discrepancies 
is  the  same  as  for  the  irradiance  reflectance,  namely  gelbstoff  fluorescence. 

4.  Conclusions 

New  parameterizations  of  the  irradiance  reflectance  and  the  remote  sensing  reflectance  in  deep 
and  shallow  waters  were  developed  using  only  the  inherent  optical  properties  of  the  water,  the 
viewing  and  solar  zenith  angle,  and  the  surface  wind  speed.  Additionally,  a  new  parameteriza¬ 
tion  for  the  upward  diffuse  attenuation  coefficient  was  developed.  The  new  model  separates  the 
dependences  on  inherent  optical  properties,  wind  speed,  viewing,  and  solar  zenith  angle.  Thus, 
their  influences  can  be  analyzed  very  easily.  The  irradiance  reflectance  and  remote  sensing  re¬ 
flectance  can  be  calculated  much  faster  using  the  analytical  equations  than  with  Hydrolight 
or  Monte  Carlo  methods.  The  estimations  of  the  irradiance  reflectance  and  remote  sensing  re¬ 
flectance  agree  significantly  better  with  the  simulations  of  Hydrolight  than  estimations  with 
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Fig.  14.  Remote  sensing  reflectance  of  shallow  water  for  the  spectral  range  from  400  to 
750  nm  for  three  different  cases  and  for  a  subsurface  viewing  angle  of  Gy  —  7°  on  the  top 
and  By  =  27°  on  the  bottom.  The  left  part  shows  the  comparison  of  simulated  (dotted  lines) 
and  estimated  values  (solid  lines)  and  the  right  side  the  relative  errors.  The  numbers  refer 
to  the  following  situations:  (1)  sediment  at  zb  —  5  m,  Cp  =  10.8  pg/1,  Cx  =  50.0  mg/1, 
ay(440nm)  =  0.2  m-1,  0V  =  27°,  u  -  1  m/s;  (2)  macrophytes  at  zb  ~  6  m,  Cp  —  2.5  fig/\, 
Cx  =  7.0  mg/1,  ay(440nm)  =  0.3  m-1,  0S  =  33°,  u  =  0  m/s;  (3)  sediment  at  zb  —  5  m, 
Cp  =  1.0  /ig/1,  Cx  =  1.0  mg/1,  ay(440nm)  =  0.05  m~l,  0S  ~  27°,  u  =  1  m/s. 


existing  equations.  The  mean  error  is  about  2-3%.  A  maximum  error  of  about  15%  occurs  at 
a  wavelength  of  685  nm  owing  to  the  fluorescence  of  chlorophyll  which  is  not  included  in  the 
system  of  equations  presented  here.  The  spectral  shape  of  the  simulations  with  Hydrolight  fits 
very  well  with  the  new  parameterizations.  The  relative  error  at  a  given  wavelength  is  below  5% 
for  the  irradiance  reflectance  and  below  10%  for  the  remote  sensing  reflectance  from  400  to 
750  nm.  Main  error  sources  are  the  fluorescence  of  chlorophyll  and  gelbstoff. 

Seasonal  changes  of  the  specific  optical  properties  of  the  water  constituents  of  Lake  Con¬ 
stance  were  investigated  previously.  The  variability  of  the  specific  absorption  of  phytoplankton 
was  analyzed  by  [19].  He  has  separated  five  different  algae  classes,  which  allow  to  model 
changes  of  the  optical  properties  of  phytoplankton  in  deep  water  from  1990  through  1992. 
The  impact  on  reflectance  spectra  is  small  compared  to  concentration  changes.  The  specific 
backscattering  of  suspended  matter  was  determined  by  [8].  He  observed  an  accuracy  of  about 
25%  for  the  estimated  concentration  of  suspended  matter  from  airborne  remote  sensing  data 
compared  to  in  situ  measurements,  indicating  a  low  variability  of  the  specific  backscattering 
coefficient  of  suspended  matter.  The  gelbstoff  exponent  S  varies  about  8%  after  [22]. 

With  the  new  parameterizations  a  set  of  equations  was  found  for  case-2  waters,  like  Lake 
Constance.  This  improves  upon  existing  equations  for  determining  the  concentration  of  the  wa¬ 
ter  constituents,  bottom  depth,  and  bottom  types  using  by  inversion  techniques.  The  analytical 
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equations  provide  a  fast  method  to  process  a  large  set  of  remote  sensing  data  from  hyperspec- 
tral  airborne  and  spacebome  sensors.  The  next  step  is  to  implement  analytical  equations  of  the 
fluorescence  of  chlorophyll  and  gelbstofif  and  to  test  the  equations  using  an  independent  dataset 
from  Lake  Constance  and  other  locations.  Inclusion  of  surface  effects  and  bidirectional  bottom 
effects  are  also  planned. 

Acknowledgments 

This  work  is  part  of  the  Special  Collaborative  Program  SFB  454  “Lake  Constance  littoral” 
funded  by  the  German  Research  Foundation  DFG.  Author  C.D.  Mobley  was  supported  by  the 
Environmental  Optics  Program  of  the  U.  S.  Office  of  Naval  Research.  Thanks  to  Peter  Gege 
and  Thomas  Heege  for  very  constructive  comments  and  suggestions. 


#3022 -$15.00  US 
(C)  2003  OSA 


Received  September  17,  2003;  Revised  October  23, 2003 
3  November  2003  /  Vol.  1 1,  No.  22  /  OPTICS  EXPRESS  2890 


Reference:  Biol  Bull.  207:  1-16.  (August  2004) 
©  2004  Marine  Biological  Laboratory 


Appendix  C 


Propagation  and  Perception  of  Bioluminescence: 
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Abstract.  Many  deep-sea  species,  particularly  crusta¬ 
ceans,  cephalopods,  and  fish,  use  photophores  to  illuminate 
their  ventral  surfaces  and  thus  disguise  their  silhouettes 
from  predators  viewing  them  from  below.  This  strategy  has 
several  potential  limitations,  two  of  which  are  examined 
here.  First,  a  predator  with  acute  vision  may  be  able  to 
detect  the  individual  photophores  on  the  ventral  surface. 
Second,  a  predator  may  be  able  to  detect  any  mismatch 
between  the  spectrum  of  the  bioluminescence  and  that  of  the 
background  light.  The  first  limitation  was  examined  by 
modeling  the  perceived  images  of  the  counterillumination 
of  the  squid  Abralia  veranyi  and  the  myctophid  fish  Cera- 
toscopelus  maderensis  as  a  function  of  the  distance  and 
visual  acuity  of  the  viewer.  The  second  limitation  was 
addressed  by  measuring  downwelling  irradiance  under 
moonlight  and  starlight  and  then  modeling  underwater  spec¬ 
tra.  Four  water  types  were  examined:  coastal  water  at  a 
depth  of  5  m  and  oceanic  water  at  5,  210,  and  800  m.  The 
appearance  of  the  counterillumination  was  more  affected  by 
the  visual  acuity  of  the  viewer  than  by  the  clarity  of  the 
water,  even  at  relatively  large  distances.  Species  with  high 
visual  acuity  (0.11°  resolution)  were  able  to  distinguish  the 
individual  photophores  of  some  counted lluminating  signals 
at  distances  of  several  meters,  thus  breaking  the  camouflage. 
Depth  and  the  presence  or  absence  of  moonlight  strongly 
affected  the  spectrum  of  the  background  light,  particularly 
near  the  surface.  The  increased  variability  near  the  surface 
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was  partially  offset  by  the  higher  contrast  attenuation  at 
shallow  depths,  which  reduced  the  sighting  distance  of 
mismatches.  This  research  has  implications  for  the  study  of 
spatial  resolution,  contrast  sensitivity,  and  color  discrimina¬ 
tion  in  deep-sea  visual  systems. 

Introduction 

Counterillumination  is  a  common  form  of  crypsis  in  the 
open  ocean  (Latz,  1995;  Harper  and  Case,  1999;  Widder, 
1999).  Its  prevalence  is  due  to  the  fact  that,  because  the 
downwelling  light  is  orders  of  magnitude  brighter  than  the 
upwelling  light,  even  an  animal  with  white  ventral  colora¬ 
tion  appears  as  a  black  silhouette  when  viewed  from  below 
(Johnsen,  2002).  This  is  particularly  disadvantageous  be¬ 
cause  an  object  is  detectable  at  a  far  greater  distance  when 
viewed  from  below  than  when  viewed  from  any  other  angle 
(Mertens,  1970;  Johnsen,  2002).  Aside  from  extremely 
transparent  tissue,  which  is  not  easy  to  achieve  in  larger 
species  with  complex  tissues,  the  way  to  overcome  this 
disadvantage  is  for  the  ventral  surface  to  emit  light  that 
matches  the  downwelling  light  in  intensity,  spectrum,  and 
angular  distribution.  Indeed,  this  solution  is  nearly  ubiqui¬ 
tous  in  nontransparent  mesopelagic  species,  particularly  in 
crustaceans,  fish,  and  squid  (Young  and  Roper,  1976;  Her¬ 
ring,  1977,  1985;  Widder,  1999). 

Counterilluminating  species  have  evolved  complex  strat¬ 
egies  to  match  the  intensity,  spectrum,  and  angular  distri¬ 
bution  of  the  downwelling  light  (Denton  et  al,  1972;  Young 
and  Mencher,  1980;  Herring,  1983;  Widder,  1999).  One 
aspect  that  is  poorly  understood,  however,  is  the  spatial 
distribution  of  the  photophores  (Young  and  Roper,  1976). 
While  some  species  ( e.g the  cookie  cutter  shark  Isistius 
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brasiliensis)  have  many  small  photophores  that  evenly  illu¬ 
minate  the  ventral  surface,  most  have  a  smaller  number  of 
isolated  photophores  that  produce  uneven  illumination  ( e.g ., 
Fig.  2d).  Thus,  even  if  the  photophores  match  the  spectrum 
and  intensity  of  the  downwelling  light  perfectly,  the  coun¬ 
terilluminator  will  be  visible  when  viewed  at  a  distance  that 
allows  these  individual  sources  to  be  discerned.  To  inves¬ 
tigate  this  problem,  the  effects  of  the  intervening  water  and 
the  viewer’s  visual  acuity  on  the  perceived  image  of  the 
counterillumination  must  be  understood. 

This  study  examines  the  effects  of  underwater  light  scat¬ 
tering  and  visual  acuity  on  the  perceived  images  of  coun¬ 
terillumination  signals.  The  effects  are  modeled  with  Monte 
Carlo  methods  and  image  transfer  theory,  using  data  col¬ 
lected  from  water  types  ranging  from  shallow  coastal  water 
to  the  deep  mesopelagic  zone  (800  m).  Three  visual  sys¬ 
tems,  with  high,  medium,  and  low  acuity,  are  also  exam¬ 
ined.  The  goal  is  to  determine  under  which  conditions 
counterilluminators  are  still  visible  and  what  implications 
this  has  for  both  camouflage  and  visual  detection  under 
low-light  conditions. 

Materials  and  Methods 

General  principles  of  image  transfer 

The  perceived  image  of  a  counterilluminating  animal 
viewed  from  a  distance  is  affected  by  three  factors:  absorp¬ 
tion  and  scattering  by  the  water  and  the  acuity  of  the 


viewer’s  eye.  The  water  and  associated  particulates  poten¬ 
tially  dim  and  blur  the  image,  and  the  acuity  of  the  eye 
determines  the  resolution  of  the  perceived  image. 

The  effect  of  the  first  factor  is  generally  modeled  in  the 
following  way.  First,  the  optical  effects  of  the  water  on  the 
image  of  a  point  source  are  calculated.  The  image  of  a  point 
source  is  known  as  the  point  spread  function  (PSF)  (Mertens 
and  Replogle,  1977).  The  point  source  is  then  convolved 
with  a  given  image  to  determine  the  appearance  of  the 
image  after  it  passes  through  the  water.  In  a  convolution, 
each  point  in  the  image  is  replaced  by  its  product  with  the 
point  spread  function  (Fig.  1).  Fortunately,  this  computa¬ 
tionally  expensive  procedure  can  be  streamlined  using  the 
convolution  theorem,  which  states  that  for  any  two  images 
/j  and  /2,  the  convolution  of  I{  with  12  is  equal  to  the  inverse 
Fourier  transform  of  the  product  of  the  Fourier  transforms 
of  the  two  images:  that  is, 

1\ *h  =  ^_1[^(/,)  *  ^(72)],  (Fig.  1)  (Equation  1) 

where  *  denotes  convolution,  and  SF(7)  and  8F-1(7)  are  the 
Fourier  and  inverse  Fourier  transforms  of  an  image  7  (Good¬ 
man,  1996).  Let  be  the  image  of  the  counterillumination, 
and  I2  be  the  point  spread  function.  Substituting  into  equa¬ 
tion  (1)  gives 

image*PSF  =  ^-l[^(image)  •  ^(PSF)].  (Equation  2) 
The  Fourier  transform  of  the  point  spread  function  is  gen- 
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Figure  1.  The  convolution  of  image  1  and  image  2  (denoted  by  the  "**’  operator)  can  be  calculated  by 
multiplying  the  Fourier  transforms  of  the  two  images  and  then  calculating  the  inverse  Fourier  transform  of  the 
product. 
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erally  referred  to  as  the  optical  transfer  function  (OTF).  Due 
to  the  convolution  theorem,  the  OTF  of  a  whole  system  is 
simply  the  product  of  the  OTFs  of  the  various  components 
in  the  system  (Goodman,  1996).  Thus,  for  this  study 

imageBnil  =  ^'‘[^(image.^J  •  OTFwater  •  OTFeyJ, 

(Equation  3) 

where  imagelinal  is  the  perceived  image  and  OTFwater  and 
OTFcyc  are  the  optical  transfer  functions  of  the  water  and 
eye  respectively.  A  final,  convenient  implication  of  the 
convolution  theorem  is  that  the  OTF  of  x  meters  of  optically 
homogeneous  water  is  equal  to  the  OTF  of  1  meter  of  water 
to  the  Xth  power.  Thus,  one  need  only  calculate  the  PSF  for 
one  distance.  This  property,  known  as  the  linearity  assump¬ 
tion,  does  not  hold  in  extreme  cases  ( e.g .,  very  large  jc),  but 
is  appropriate  for  the  situations  examined  in  this  study 
(Jaffe,  1992).  The  equation  for  underwater  image  transfer  is 
then 

imageHnal(x)  =  9;“l[^(imagem,IiJll)  •  (OTF,)*  •  OTFtye], 

(Equation  4) 

where  imagertnal(x)  is  the  perceived  image  viewed  from  a 
distance  of  x  meters,  and  OTF!  is  the  optical  transfer  func¬ 
tion  of  the  water  over  a  distance  of  1  meter. 

Although  Eq.  (4)  correctly  describes  the  propagation  of  a 
two-dimensional  image,  it  requires  modification  when  used 
in  the  context  of  counterillumination,  because  the  back¬ 
ground  radiance  is  affected  by  the  entire  three-dimensional 
light  field  and  changes  as  the  viewer  moves  down  and  away 
from  its  target.  From  Mertens  (1970),  the  degradation  of 
contrast  of  a  large  image  underwater  ( i.e .,  the  OTF  at  zero 
spatial  frequency)  is 


Cr 

OTF(O)  =  -£■  =  e~(c~KOx, 
Co 


(Equation  5) 


where  Cx  and  C0  are  contrast  at  x  and  0  meters  viewing 
distance,  c  is  the  beam  attenuation  coefficient,  and  KL  is  the 
attenuation  coefficient  of  the  background  radiance.  In  the 
case  of  upward  viewing,  KL  equals  Kl*  the  attenuation 
coefficient  of  direct  downward  radiance  (Johnsen,  2002). 

The  correct  OTF  for  objects  being  viewed  from  below  is 
obtained  by  normalizing  the  original  OTF  so  that  OTF(O) 
equals  e~(c  ~  Kut)x  (Mertens,  1970).  Thus  the  final  equation 
for  the  propagation  of  images  viewed  from  below  is 


image^M 


=  9?-' 


^(image^J  • 


OTF,  V 
OTF,(0)j  '°TF< 


(Equation  6) 


The  OTF  is  a  complex-valued  function  and  difficult  to 
interpret.  Therefore,  its  magnitude,  known  as  the  modula¬ 


tion  transfer  function  (MTF),  is  often  also  calculated.  The 
MTF  is  quite  useful  because  it  gives  the  fraction  of  remain¬ 
ing  image  contrast  as  a  function  of  spatial  frequency.  For 
example,  MTF(4)  =  0.5  implies  that  50%  of  the  original 
image  contrast  remains  for  detail  that  has  a  spatial  fre¬ 
quency  of  4  cycles  per  degree. 


Images  examined 

Images  of  the  ventral  bioluminescence  of  two  counter- 
illuminating  species  were  used:  (1)  the  enoploteuthid  squid 
Abralia  veranyi  (Riippell,  1844)  (eye-flash  squid),  and  (2) 
the  myctophid  fish  Ceratoscopelus  maderensis  (Gunther, 
1864)  (horned  lanternfish)  (Fig.  2A,  B).  The  two  were 
chosen  to  provide  a  range  of  photophore  spacing.  Counter¬ 
illumination  in  A.  veranyi  is  finely  detailed;  that  of  C. 
maderensis  is  relatively  coarse  (Fig.  2C,  D).  A.  veranyi  was 
collected  at  depth,  using  the  Johnson-Sea-Link  research 
submersible  fitted  with  1 1 -liter  acrylic  plastic  cylinders  with 
hydraulically  activated,  sliding  lids.  C.  maderensis  was  col¬ 
lected  at  night,  using  an  opening/closing  Tucker  trawl 
(4.3-m2  opening,  lA  inch  knotless  nylon  mesh)  fitted  with  a 
thermally  insulated  collecting  container.  Specimens  were 
manually  stimulated  to  bioluminesce,  and  then  were  re¬ 
corded  with  a  Dage  ISIT  image-intensified  video  camera  (A. 
veranyi)  or  Intevac  GenllSys  image  intensifier  system  and 
CCD  video  camera  (C.  maderensis).  Images  that  show  how 
the  counterilluminating  animals  appear  from  below  (Fig. 
2E,  F)  were  created  by  combining  the  bioluminescence 
images  with  silhouettes  of  the  animals  obtained  from  nor¬ 
mal  illumination  photographs  (taken  immediately  after  the 
intensified  images).  Non-illuminating  portions  of  the  ani¬ 
mals  are  shown  as  black  because  this  is  how  they  appear 
against  the  downwelling  light  (Johnsen,  2002).  The  natural 
posture  of  A.  veranyi  is  unknown.  Although  observers  in 
submersibles  generally  find  mesopelagic  squid  with  their 
fins  folded  and  their  arms  and  tentacles  placed  over  their 
heads  (Vecchione  and  Roper,  1991;  Fig.  2A),  this  is  may  be 
a  response  to  the  perceived  threat  from  the  submersible.  In 
the  silhouette  chosen,  the  fins  and  appendages  were  ex¬ 
tended  to  examine  their  effect  on  visibility. 

The  backgrounds  were  set  to  a  brightness  equal  to  the 
average  brightness  of  the  counterilluminating  animal  (minus 
the  fins,  arms,  and  tentacles  in  the  case  of  the  squid). 
Because  these  relative  values  allow  the  animal  to  blend  with 
the  background  most  easily,  it  is  assumed  that  they  approx¬ 
imately  match  what  would  be  observed  in  the  field.  The 
backgrounds  for  the  C.  maderensis  images  are  darker  be¬ 
cause  the  average  brightness  of  the  animal  is  darker  (due  to 
the  wider  spacing  of  the  photophores).  Note  that  these 
figures  show  relative  brightnesses,  chosen  to  maximize  vis¬ 
ibility  on  the  printed  page.  The  absolute  brightnesses  are  of 
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Figure  2.  (A)  The  eye-flash  squid  Ahralia  veranyi.  (B)  The  homed  lantemfish  Ceratoscopelus  maderensis. 

(C)  Counterillumination  of  A.  veranyi.  (D)  Counterillumination  of  C.  maderensis.  (E)  Counterillummation  of  A. 
veranyi  as  viewed  from  below  against  the  downwelling  light.  (F)  Counterillumination  of  C.  maderensis  as 
viewed  from  below  against  the  downwelling  light.  Scale  bar  is  5  cm.  Background  in  (E)  and  (F)  is  set  to  the 
average  brightness  of  the  counterilluminating  animal.  Panel  B  courtesy  of  Marine  Biological  Laboratory  Digital 
Archive,  Flescher  Fish  Collection. 


course  much  dimmer  (far  beyond  the  reach  of  printed  paper) 
and  can  only  be  seen  by  the  dark-adapted  eye. 

The  intensified  images  are  not  perfect  representations  of 
the  actual  counterillumination.  The  resolution  of  the  images 
is  low,  and  the  photophore  signals  are  slightly  expanded  due 
to  “blooming”  of  the  image  at  the  detector  array.  In  addi¬ 
tion,  although  counterillumination  is  more  stable  than  other 
bioluminescent  signals,  the  images  are  static  representations 
of  potentially  variable  light  emission.  Indeed,  a  subset  of  the 
ventral  photophores  in  A.  veranyi  was  not  lit  in  the  studied 
image  (Herring  et  al. ,  1992).  This  relatively  low  number  of 
small  photophores  most  likely  would  not  change  a  spatial 
distribution  that  is  already  quite  uniform.  However,  they 
may  play  a  role  in  spectral  changes.  In  C.  maderensis,  all  the 
ventral  photophores  were  emitting  during  the  image  expo¬ 
sure. 


Calculation  of  point  spread  functions 
and  attenuation  coefficients 

The  PSFs  in  this  study  were  determined  using  Monte 
Carlo  software  (BSFPSF  ver.  1.1.,  developed  by  CDM). 
Five  million  simulated  photons  were  tracked  from  an  iso¬ 
tropic  point  source  (of  unit  power)  to  their  point  of  inter¬ 
section  with  a  sphere  of  radius  1  m.  Although  a  PSF  is 
defined  as  the  image  of  a  cosine  point  source,  the  use  of  an 
isotropic  point  source  achieves  the  same  result  because 
scattering  in  natural  waters  is  primarily  in  a  forward  direc¬ 
tion  (Mertens  and  Replogle,  1977;  confirmed  by  prelimi¬ 
nary  calculations).  Due  to  the  symmetry  of  an  isotropic 
point  source,  calculations  could  be  completed  in  far  less 
time  than  if  a  cosine  point  source  were  used. 

The  radiance  distribution  of  the  simulated  photons  at  the 
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intersection  with  the  1-m  sphere  is  the  PSF.  The  three 
factors  (besides  distance)  affecting  the  PSF  are  (1)  the 
absorption  coefficient  a,  (2)  the  scattering  coefficient  b,  and 
(3)  the  phase  function  y.  The  first  and  second  factors  specify 
how  often  a  photon  is  absorbed  or  scattered  by  the  water  and 
associated  particulates.  The  third  factor  specifies  the  angular 
distribution  of  the  scattered  light.  Absorption  and  attenua¬ 
tion  coefficients  were  obtained  for  four  water  types:  (1) 
coastal  water  at  5-m  depth,  (2)  oceanic  water  at  5-m  depth, 
(3)  oceanic  water  at  210-m  depth,  and  (4)  oceanic  water  at 
800-m  depth  (Table  1).  Absorption  and  scattering  coeffi¬ 
cients  for  coastal  water  were  obtained  by  Dr.  Heidi  Sosik 
(Woods  Hole  Oceanographic  Institution,  Woods  Hole,  MA) 
using  a  dual-path,  multiband  absorption/attenuation  meter 
(ac-9,  WETLabs  Inc.)  at  a  site  80  km  from  the  coast  of 
Portsmouth,  New  Hampshire  (42°47'N  70°05'W,  1106  lo¬ 
cal  time,  30  June  2000)  (see  Johnsen  and  Sosik,  2003,  for 
details).  Optical  coefficients  in  oceanic  water  (Jerlov  type  I) 
at  5  m  and  210  m  were  obtained  by  Drs.  Andrew  Barnard, 
Scott  Pegau,  and  Ronald  Zaneveld  (College  of  Oceanic  and 
Atmospheric  Sciences,  Oregon  State  University,  Corvallis, 
Oregon)  using  the  same  equipment  in  the  equatorial  Pacific 
(1005  local  time,  30  April  1996;  0°0'N  177°21'W).  Optical 
coefficients  in  oceanic  water  at  800-m  depth  were  obtained 
from  Capone  et  al.  (2002).  In  all  cases,  absorption  and  beam 
attenuation  coefficients  were  measured  at  412,  440,  488, 
510,  532,  555,  650,  and  676  nm.  Although  point  spread 
functions  were  calculated  for  all  eight  wavelengths,  for 
clarity  only  those  at  412,  488,  555,  and  650  nm  are  analyzed 
and  discussed  in  this  study.  While  the  5-m  coastal  measure¬ 
ment  is  somewhat  specific  to  measurement  site,  the  three 
oceanic  measurements  are  typical  of  most  oceanic  waters, 
particular  those  at  210  and  800  m. 

Because  the  ac-9  absorption-attenuation  meter  has  detec¬ 
tors  of  finite  size,  light  scattered  over  small  angles  was 
collected  by  the  detector  and  incorrectly  interpreted  as  un¬ 
scattered.  Thus,  scattering  was  underestimated  by  a  small 
amount.  If  one  assumes  that  the  scattering  matches  Pet- 
zold’s  phase  function,  then  the  scattering  coefficient  is  un¬ 
derestimated  by  approximately  20%.  Again,  preliminary 
results  showed  that  this  had  negligible  effect  on  the  blurring 
of  the  image,  though  it  would  have  resulted  in  slightly 


greater  attenuation  of  the  contrast  of  the  whole  image.  The 
ac-9  meter  also  does  not  measure  certain  factors  that  may 
influence  image  propagation,  such  as  marine  snow  and 
refractive  index  inhomogeneities.  The  large  particles  of 
marine  snow  will  limit  the  long-range  visibility  of  small 
objects  by  direct  occlusion,  and  refractive  index  inhomoge¬ 
neities  may  slightly  increase  scattering  at  very  small  angles 
(below  the  resolution  limit  of  the  visual  systems  examined) 
(Bogucki  et  al.,  1998). 

The  phase  function  y  was  chosen  to  be  the  commonly 
used  “average  particle”  function  (Mobley  et  al.,  1993)  based 
on  measurements  by  Petzold  (1977).  In  productive  coastal 
waters,  most  of  the  light  is  scattered  by  living  phytoplank¬ 
ton,  which  have  a  backscatter  fraction  of  0.01  or  less  ( e.g ., 
Ulloa  et  al.,  1994).  However,  in  clear  oceanic  water,  iso¬ 
tropic  scattering  by  the  water  itself  is  a  significant  fraction 
of  the  total  scattering,  and  the  total  backscatter  fraction  can 
be  as  large  as  0.04  (Mobley,  1994).  We  chose  Petzold’s 
average  particle  phase  function  (Mobley  et  al.,  1993),  which 
has  a  backscatter  fraction  of  0.018,  about  midway  between 
the  two  extremes.  Preliminary  results  showed  that,  because 
scattered  light  was  extremely  dim  compared  to  the  unscat¬ 
tered  direct  beam,  the  choice  of  phase  function  made  no 
notable  difference. 

PSF  values  were  calculated  up  to  10°,  at  0.05°  intervals. 
Although  the  PSF  from  0  to  1  °  was  calculated  using  Monte 
Carlo  methods,  computational  limits  (due  to  the  small  size 
of  the  angular  bins  receiving  scattered  photons)  prevented 
accurate  calculations  at  substantially  higher  angles  for  the 
given  number  of  initial  photons.  Therefore,  the  PSF  from  1° 
to  10°  was  estimated  by  fitting  the  PSF  from  0.45°  to  1°  to 
a  power  function  and  then  extrapolating  by  0.05°  incre¬ 
ments  up  to  an  angle  of  10°  (see  Voss,  1991,  for  justifica¬ 
tion). 

The  optical  transfer  functions  of  the  eyes  were  modeled 
as  the  Gaussian  curve: 

OTF(v)  =  e~356(R^3  (Equation  7) 

where  v  is  the  spatial  frequency  (in  cycles/degree)  and  R  is 
the  spatial  resolution  (Warrant,  1999).  This  function,  often 
used  to  model  the  OTF  of  visual  systems,  results  in  a  barely 


Table  1 

Absorption  and  scattering  coefficients  (a  and  b  respectively)  used  in  the  Monte  Carlo  calculation  of  point  spread  functions 


Coastal  water  Oceanic  water  Oceanic  water  Oceanic  water 

at  5-m  depth  at  5-m  depth  at  210-m  depth  at  800-m  depth 

Wavelength 

(nm)  a  b  a  b  a  h  a  h 


412  0.29  0.26  0.035  0.11  0.060  0.018  0.027  0.020 

488  0.15  0.21  0.038  0.098  0.035  0.013  0.027  0.014 

555  0.11  0.19  0.073  0.091  0.077  0.0094  0.072  0.0060 

650  0.37  0.15  0.36  0.085  0.36  0.014  0.36  0.014 


6 


S.  JOHNSEN  ET  AL. 


detectable  contrast  of  2.8%  (=  e~3  56)  at  the  spatial  resolu¬ 
tion  of  the  eye.  The  spatial  resolutions  of  three  mesopelagic 
fish  were  chosen  to  span  a  wide  range  of  visual  acuity:  (1) 
R  =  0.11°  (the  “lovely  hatchetfish’’  Argyropelecus  aculea- 
tus ),  (2)  R  =  0.23°  (the  spookfish  Opisthoproctus  soleatus ), 
and  (3)  R  =  0.50°  (the  myctophid  fish  Lampanyctus  festi- 
vus)  (Collin  et  al.,  1997;  Wagner  et  al,  1998).  A.  aculeatus 
and  O.  soleatus  both  have  up  ward- viewing  tubular  eyes;  L. 
festivus  has  lateral-viewing  eyes  and  so  probably  does  not 
search  for  overhead,  counterilluminating  prey. 

The  acuity  of  these  species  was  measured  from  the  den¬ 
sity  of  their  retinal  ganglion  cells  (which  accounts  for  spa¬ 
tial  summation).  Because  these  counts  also  include  dis¬ 
placed  ganglion  cells,  they  may  slightly  overestimate  acuity. 
The  predicted  acuity  also  assumes  a  well-focused  image,  but 
this  is  generally  the  case  for  the  foveal  regions  of  deep-sea 
eyes  (Warrant  and  Locket,  2004).  Increasing  spatial  sum¬ 
mation  will  also  lower  the  acuity.  Finally,  it  is  important  to 
note  that  these  spatial  resolutions  do  not  include  potential 
blurring  of  a  moving  image  due  to  large  temporal  summa¬ 
tion.  Since  long  temporal  summation  times  are  common  at 
depth  (Frank,  1999)  and  animals  do  drift  relative  to  one 
another  in  the  water,  the  actual  spatial  resolution  in  certain 
situations  may  be  less  than  that  predicted  by  retinal  mor¬ 
phology. 

The  minimum  contrast  threshold  is  the  smallest  percent¬ 
age  variation  in  radiance  that  can  be  detected.  This  value  for 
fish  is  approximately  l%-2%  in  bright  light,  but  it  rises  as 
depth  increases  (Douglas  and  Hawryshyn,  1990).  Though 
few  direct  measurements  have  been  made,  the  threshold  at 
mesopelagic  light  levels  appears  to  range  from  about  25%  to 
50%  ( e.g .,  threshold  for  the  Atlantic  cod  Gadus  morhua  at 
650-m  depth  is  approximately  50%  (Anthony,  1981)).  We 
therefore  set  the  minimum  contrast  threshold  at  33%,  while 
accepting  that  depth,  water  clarity,  and  special  visual  adap¬ 
tations  make  the  actual  threshold  highly  variable. 

The  attenuation  coefficients  of  direct  downward  radiance 
Kul  were  calculated  by  modeling  the  underwater  radiance 
distribution  using  radiative  transfer  software  (Hydrolight  4.2, 
Sequoia  Scientific  Inc.,  Bellevue,  WA,  www.hydrolight.info). 
The  inherent  optical  properties  required  by  the  software  were 
obtained  from  measured  vertical  profiles  of  chlorophyll  con¬ 
centration  and  absorption  and  scattering  coefficients  from  the 
four  water  types  examined  (see  Johnsen,  2002;  and  Johnsen 
and  Sosik,  2003,  for  details).  The  sun  was  assumed  to  be  at  the 
zenith  on  a  clear  day  with  no  wind. 

Measurement  of  moonlight  and  starlight  spectra 

Nocturnal  spectra  under  moonlight  and  starlight  were 
measured  using  a  spectrometer  with  a  highly  sensitive  pho¬ 
tomultiplier  detector  (OL-754-PMT,  Optronics  Laboratories 


Inc.,  Orlando,  FL).  Moonlight  spectra  were  measured  in  air 
on  a  barrier  island  in  Florida  during  full  moon  (moon  was  at 
its  peak  elevation).  An  integrating  sphere  was  used  to  col¬ 
lect  light  from  all  regions  of  the  sky.  Starlight  spectra  were 
measured  on  a  moonless  night  on  a  completely  darkened 
ship  in  the  center  of  the  Gulf  Stream  (latitude  ~27°N)  to 
ensure  a  complete  absence  of  light  pollution.  To  minimize 
light  loss,  the  integrating  sphere  was  removed  and  the 
entrance  slit  of  the  spectrometer  (—30°  angular  field)  was 
aimed  at  the  zenith.  The  downwelling  irradiance  at  5-m 
depth  under  moonlight  and  starlight  was  calculated  using 
the  above-described  radiative  transfer  software,  with  the 
correct  skylight  spectrum  as  an  input. 

Results 

Point  spread  and  optical  transfer  functions  of  the  water 
types  and  visual  resolutions 

The  point  spread  functions  in  all  four  water  types  were 
extremely  narrow,  with  the  radiance  at  zero  degrees  2-3 
orders  of  magnitude  larger  than  the  radiance  at  higher 
angles  (at  a  distance  of  5  m)  (Fig.  3).  With  increasing  water 
clarity  and  depth,  this  effect  became  more  pronounced.  The 
wavelength  dependence  of  the  PSF  was  complex,  depend¬ 
ing  on  the  relative  numbers  of  scattering  and  absorption 
events. 

In  all  four  water  types,  the  modulation  transfer  function 
was  primarily  affected  by  the  visual  resolution  of  the  view¬ 
er’s  eye  (Fig.  4).  However,  the  MTFs  in  near-surface  waters 
decreased  at  higher  spatial  frequencies  (independently  of 
the  decrease  due  to  visual  acuity  limitations),  indicating 
some  blurring  by  the  water  (Fig.  4A,  B).  The  MTFs  within 
a  given  water  type  had  similar  shapes  and  differed  primarily 
in  magnitude  (set  by  MTF(0)  =  e~5(c  ~  Ku)).  This  magni¬ 
tude  had  a  complicated  wavelength  dependence,  being  pro¬ 
portional  to  wavelength  in  near- surface  waters  and  inversely 
proportional  to  wavelength  in  deep  waters. 

Perceived  images 

The  perceived  images  were  dramatically  affected  by  the 
visual  resolution  of  the  viewer  and,  to  a  lesser  extent,  by 
scattering  and  absorption  by  the  water  (Figs.  5,  6).  When 
viewed  at  0.11°  resolution,  Ceratoscopelus  maderensis  and 
Abralia  veranyi  had  a  contrast  greater  than  33%  at  distances 
up  to  4  to  8  m  (though  the  visibility  of  the  latter  was 
primarily  due  to  the  unlit  fins  and  appendages).  However, 
when  viewed  at  0.5°  resolution,  the  contrast  of  the  counter¬ 
illumination  was  greater  than  33%  only  up  to  a  distance  of 
1  to  2  m.  The  individual  photophores  of  C.  maderensis  were 
distinguishable  up  to  2  m  at  0.1 1°  resolution,  and  the  gen¬ 
eral  pattern  of  photophores  was  distinguishable  up  to  2  m  at 
lower  resolutions.  The  general  pattern  of  the  photophores  of 
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Figure  3.  Radiance  vs.  angle  for  a  point  source  viewed  from  a  distance  of  5  in  (point  spread  function).  The 
radiance  is  normalized  by  the  radiance  of  a  point  source  viewed  at  5  m  in  a  medium  that  does  not  scatter  or 
absorb  light.  (A)  Coastal  water  at  5-m  depth  (B)  Oceanic  water  at  5-m  depth.  (C)  Oceanic  water  at  210-m  depth. 
(D)  Oceanic  water  at  800-m  depth.  The  normalized  radiances  at  zero  degrees  are  given  numerically  rather  than 
graphically  because  they  are  far  higher  than  the  other  values. 


A.  veranyi  was  evident  at  1  m  at  0.1 1°  resolution,  but  not  at 
lower  resolutions. 

The  counterillumination  of  both  species  is  visible  at 
roughly  twice  the  distance  in  the  clearest  conditions  studied 
(at  488  nm  in  oceanic  water  at  800-m  depth)  as  it  is  in  the 
most  turbid  conditions  (at  412  nm  in  coastal  water  at  5-m 
depth)  (Fig.  5A  vs.  5B,  Fig.  6A  vs.  6B).  The  difference  was 
entirely  due  to  the  difference  in  MTF(O)  between  the  two 
water  types  (95%  vs.  37%  at  a  distance  of  5  m)  and  not  to 
significantly  increased  blurring  of  fine  detail. 

Variation  of  background  spectra  and  wavelength 
dependence  of  contrast  attenuation 

The  background  spectra  at  shallow  depths  under  moon¬ 
light  and  starlight  differed  substantially  in  both  coastal  and 
oceanic  waters,  particularly  at  shorter  wavelengths  (Fig.  7A, 
B).  Under  starlight,  the  spectrum  narrowed  and  the  peak 
wavelength  was  long- shifted  (by  40  to  80  nm  depending  on 


the  water  type  and  what  is  considered  the  true  peak).  The 
background  spectra  were  also  affected  substantially  by 
depth,  even  at  mesopelagic  depths.  As  the  depth  increased 
from  200  to  800  m,  the  spectra  of  the  downwelling  irradi- 
ance  narrowed  slightly  and  the  peak  wavelength  shifted 
from  490  nm  to  480  nm  (Fig.  7C). 

General  contrast  attenuation  was  relatively  rapid  and 
wavelength-independent  at  5-m  depth  in  both  coastal  and 
oceanic  waters,  with  sighting  distances  (proportional  to 
l/c  —  Ku)  only  5%-20%  of  those  in  deeper  waters  (Fig. 
7D).  At  greater  depths,  sighting  distance  was  highly  depen¬ 
dent  on  wavelength.  At  these  depths,  sighting  distance  in¬ 
creased  with  wavelength,  until  it  reached  a  peak  at  a  wave¬ 
length  about  30  nm  longer  than  that  of  the  peak  wavelength 
of  downwelling  irradiance.  After  this  peak,  the  sighting 
distance  decreased  rapidly  with  wavelength.  For  wave¬ 
lengths  greater  than  600  nm,  the  sighting  distances  at  depth 
were  less  than  those  near  the  surface. 
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Figure  4.  Contrast  as  a  function  of  spatial  frequency  for  an  object  viewed  from  a  distance  of  5  m 
(modulauon  transfer  function,  MTF).  The  contrast  is  normalized  by  the  contrast  at  a  distance  of  0  m.  (A)  Coastal 
water  at  5-m  depth.  (B)  Oceanic  water  at  5-m  depth.  (C)  Oceanic  water  at  210-m  depth.  (D)  Oceanic  water  at 
800-m  depth.  The  MTF  is  shown  for  two  visual  systems,  one  with  0.1 1°  resolution  and  one  with  0.5°  resolution. 

At  a  spatial  frequency  of  approximately  0  4  cycles/deg,  the  data  split,  with  the  lower  trace  denoting  0.5° 
resolution.  The  two  gray  lines  in  (A)  denote  the  MTF  for  the  eyes  alone.  Because  the  MTF  at  0  cycles/deg  is 
important,  the  graphs  include  this  point  despite  bemg  logarithmic. 


Discussion 

This  study  shows  that  a  counterilluminating  individual 
faces  a  number  of  difficulties.  First,  an  acute  eye  (0.11° 
resolution)  with  moderate  contrast  sensitivity  (33%)  can 
detect  the  photophore  patterns  of  both  Abralia  veranyi  and 
Ceratoscopelus  maderensis  at  distances  greater  than  1  m. 
Second,  even  the  water  at  the  relatively  turbid  shallow 
coastal  site  blurred  the  counterillumination  signals  very 
little.  Although  all  four  water  types  did  lower  the  overall 
contrast  of  the  counterilluminator,  the  attenuation  rate  was 
quite  low,  particularly  at  mesopelagic  (>200  m)  depths. 
Finally,  the  spectrum  of  down  welling  background  light  var¬ 
ied  considerably  with  depth  in  the  mesopelagic  zone  and 
was  strongly  affected  by  the  source  of  nocturnal  illumina¬ 
tion  at  the  shallow  depths.  This  suggests  that  counterillu¬ 
minating  photophores  must  be  spaced  closely  together  when 
viewed  by  visually  acute  species,  and  that  matching  the 


background  spectrum  may  be  more  difficult  than  previously 
considered.  From  the  predator’s  point  of  view,  this  study 
suggests  that  high  spatial  resolution  and  color  discrimina¬ 
tion  in  the  blue-green  portion  of  the  spectrum  are  highly 
advantageous  for  detecting  counterilluminating  prey.  How¬ 
ever,  since  both  of  these  characteristics  reduce  sensitivity, 
they  also  have  a  cost  that  must  be  balanced. 

The  remainder  of  the  paper  explores  these  limitations  in 
detail.  It  is  important  to  note  that,  despite  these  limitations, 
counterillumination  dramatically  decreases  the  visibility  of 
the  individual.  The  visibility  of  A.  veranyi  at  distances 
greater  than  1  m  is  entirely  due  to  the  unlit  fins,  tentacles, 
and  arms  (Fig.  5),  which  may  be  held  above  and  against  the 
body  to  minimize  their  silhouette  (Fig.  2A).  The  visibility  of 
these  unlit  regions  at  distances  of  at  least  8  m  highlights  the 
impressive  crypsis  afforded  by  counterillumination.  In  ad¬ 
dition,  in  certain  cases  the  goal  may  not  be  complete  cryp- 
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Figure  5.  Counterillumination  of  Abralia  veranyi  viewed  from  distances  of  1,  2,  4,  and  8  m  by  animals  with 
eyes  of  0. 1 1°,  0.23°,  and  0.5°  resolution.  (A)  Countenllumination  is  viewed  at  a  wavelength  of  412  nm  in  coastal 
water  at  5-m  depth  (the  optical  conditions  that  had  the  greatest  effect  on  image  propagation).  (B)  Counterillu¬ 
mination  is  viewed  at  a  wavelength  of  488  nm  in  oceanic  water  at  800-m  depth  (the  optical  conditions  that  had 
the  least  effect  on  image  propagation).  The  percentages  indicate  the  maximum  contrast  in  each  image.  All  images 
are  scaled  in  size  for  viewing  distance,  and  the  backgrounds  are  all  set  equal.  To  see  the  absolute  brightness 
values  in  the  image,  view  the  figure  under  dim  illumination  so  that  the  printed  background  matches  the  brightness 
at  the  relevant  depth.  For  example,  to  see  what  the  counterillummation  looks  like  at  depths  of  200,  300,  and 
400  m,  view  the  figure  under  civil  twilight,  full  moonlight,  and  half-moonlight  respectively. 


sis,  but  a  bioluminescent  analog  of  disruptive  coloration. 
The  individual  photophores  may  break  up  the  silhouette  so 
that  it  appears  as  a  number  of  small  objects  rather  than  as 


one  large,  recognizable  outline.  This  tactic  is  common  and 
highly  successful  in  benthic  and  terrestrial  habitats  where 
the  background  is  complex  (Lythgoe,  1979).  Its  effective- 
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Figure  6.  Identical  to  Figure  5,  except  that  the  counterillumination  signal  is  generated  by  Cercitoscopelux 
maderensis. 


ness  in  pelagic  environments,  where  the  background  is  very 
uniform,  is  uncertain.  Finally,  the  ability  of  the  predator  to 
recognize  the  perceived  image  as  potential  prey  depends  on 
pattern  recognition,  a  higher  level  of  visual  processing  that 
is  poorly  understood  in  oceanic  species. 


Effects  of  intervening  water  on  counter  illumination 

Despite  the  authors’  initial  expectations,  the  water  had 
little  effect  on  the  appearance  of  the  counterillumination. 
This  was  due  to  several  factors.  First,  even  in  the  worst  case 
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Figure  7.  (A)  Downwelling  irradiance  (normalized  by  peak)  at  5-m  depth  in  coastal  water  under  moonlight 

and  starlight.  (B)  Same  as  (A),  but  in  oceanic  water.  (C)  Downwelling  daytime  irradiance  at  depths  of  200  and 
800  m.  (D)  Sighting  distance  of  a  counterilluminator,  normalized  by  sighting  distance  in  mesopelagic  waters  at 
480  nm.  Because  all  published  bioluminescence  spectra  are  calibrated  in  energy  units  (instead  of  quanta),  (A), 

(B)  and  (C)  are  calibrated  in  relative  energy  units. 


(412-nm  light  in  coastal  water  at  5-m  depth),  8%  of  the 
photons  were  neither  scattered  nor  absorbed  after  traveling 
5  m  ( =e~5(a  +  h)  =  e~5(°29  +  026)).  xhe  viewer  perceives 
these  uncollided  photons  as  all  arriving  from  one  point 
(which  has  a  radius  determined  by  the  spatial  resolution  of 
the  eye),  so  the  perceived  radiance  of  this  unscattered  light 
is  quite  high.  In  contrast,  the  scattered  photons,  although 
more  numerous  (92%  in  this  case),  are  scattered  over  a 
larger  angular  area  and  so  have  a  lower  radiance.  Second, 
for  a  photon  to  contribute  to  image  blurring,  it  must  be 
scattered  but  not  absorbed  by  the  water.  Because  highly 
scattering  waters  also  tend  to  be  highly  absorbing  (see  Table 
1),  many  scattered  photons  are  absorbed  before  they  can 
reach  the  eye  and  thus  cannot  contribute  to  image  blurring. 
This  is  in  sharp  contrast  to  atmospheric  scattering  ( e.g .,  due 
to  haze  or  fog),  which  can  be  high  while  absorption  is 
negligible  (Bohren  and  Huffman,  1998).  Finally,  light  scat¬ 
tering  in  water  is  strongly  biased  in  the  forward  direction.  In 
the  phase  function  used  in  this  study,  more  than  50%  of  the 
scattered  photons  are  deflected  less  than  5°  (Petzold,  1977). 


Thus,  the  halo  of  scattered  light  surrounding  the  image  of  a 
point  source  is  quite  narrow.  This  also  differs  from  the 
atmospheric  case,  where  light  is  often  scattered  over  rela¬ 
tively  large  angles  (Bohren  and  Huffman,  1998).  For  all  the 
above  reasons,  the  images  in  all  the  water  types  examined 
did  not  lose  a  substantial  amount  of  fine  detail.  There  are,  of 
course,  considerably  more  turbid  marine  waters  (very  close 
to  shore  or  to  river  plumes,  coccolithophore  blooms,  etc.). 
Counterilluminators,  however,  are  seldom  found  in  these 
locations. 

A  parameter  that  was  greatly  affected  by  the  water  was 
the  attenuation  of  the  contrast  of  the  entire  scene  (e.g., 
MTF(0)).  The  attenuation  coefficient  depends  on  the  view¬ 
ing  angle  of  the  predator  and  for  upward  viewing  is  c  — 

(see  Eq.  5).  This  coefficient  is  far  smaller  than  the  attenu¬ 
ation  coefficients  for  viewing  an  object  horizontally  or  from 
above  (c,  and  c  +  Ku  respectively),  so  objects  viewed  from 
below  are  visible  at  much  longer  distances  (Johnsen,  2002). 
This  result  derives  from  the  fact  that,  as  a  viewer  moves 
down  and  away  from  a  counterilluminator,  the  background 
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dims  almost  as  quickly  as  the  signal  does,  thus  maintaining 
the  contrast.  The  unusual  wavelength  dependence  of  the 
attenuation  of  counterillumination  occurs  because  c  and  K ^ 
vary  somewhat  independently.  At  shorter  wavelengths,  c  — 
Hu  increases  slightly  with  wavelength;  at  longer  wave¬ 
lengths  at  depth,  c  increases  rapidly  with  wavelength,  while 
Ku  remains  more  or  less  constant.  This  is  because  almost 
all  long-wavelength  light  at  depth  is  due  to  Raman  scatter¬ 
ing,  in  which  a  small  portion  of  the  blue-green  light  is 
converted  into  longer  wavelength  light  (Marshall  and 
Smith,  1990;  Johnsen,  2002).  Because  this  Raman-scattered 
light  is  produced  de  novo  from  shorter  wavelength  light,  it 
has  roughly  the  same  attenuation  coefficient  as  that  light, 
and  so  the  difference  between  c  and  K ^  can  grow  quite 
large.  But  because  the  long  wavelength  light  is  extremely 
dim,  it  may  not  be  of  visual  significance,  particularly  at 
mesopelagic  depths. 

A  curious  feature  of  this  wavelength  dependence  is  that 
the  wavelength  of  least  contrast  attenuation  is  about  30  nm 
longer  than  the  peak  wavelength  of  the  downwelling  light. 
The  lower  contrast  attenuation  at  these  wavelengths  allows 
for  a  slightly  longer  sighting  distance  (proportional  to 
1/c  —  Ku\  12.5%  longer  at  210  m;  5.5%  longer  at  800  m) 
than  at  the  peak  wavelength.  Because  the  spectral  responses 
of  most  deep-sea  visual  systems  are  relatively  flat  (Douglas 
et  al.,  1998),  this  shift  may  be  inconsequential. 

Effect  of  variation  in  background  illumination 

The  fact  that  the  spectrum  of  the  background  changes 
with  depth  has  been  examined  before  (e.g.,  Young  and 
Mencher,  1980).  This  study  confirms  that,  even  at  mesope¬ 
lagic  depths,  the  spectrum  changes  substantially  with  depth. 
While  a  10-nm  shift  in  the  peak  wavelength  appears  minor, 
it  causes  large  shifts  in  the  intensity  of  the  off-peak  light 
because  the  wavelength  distributions  are  quite  narrow.  For 
example,  if  the  peak  intensities  are  set  equal  at  100%  (as  in 
Fig.  7C),  the  difference  between  the  downwelling  irradiance 
at  depths  of  200  and  800  m  is  62%  at  500  nm  and  32%  at 
470  nm. 

A  previously  unconsidered  issue  is  the  effect  of  the 
nocturnal  illumination  source.  Many  counterilluminators 
are  vertical  migrators  and  can  be  found  near  the  surface  at 
night  (the  downwelling  irradiance  at  5-m  depth  under 
moonlight  and  starlight  equals  that  found  during  the  middle 
of  the  day  at  300  and  500  m  respectively).  The  background 
illumination  then  depends  on  whether  the  moon  is  present. 
Over  a  complete  lunar  cycle,  the  moon  is  above  the  horizon 
for  about  half  of  the  nocturnal  hours.  Because  the  moon 
reflects  all  wavelengths  more  or  less  equally  (Munz  and 
McFarland,  1977),  the  spectrum  of  the  night  sky  with  the 
moon  present  is  similar  to  the  spectrum  of  daylight  (though 
dimmer  by  about  6  orders  of  magnitude,  and  slightly  red- 
shifted  due  to  background  starlight).  When  the  moon  is  not 


present,  the  illumination  has  three  primary  components:  (1) 
starlight,  mostly  due  to  dim  red  stars  invisible  to  the  naked 
eye,  (2)  scattering  of  sunlight  by  dust  in  the  plane  of  the 
solar  system,  and  (3)  emission  spectra  from  gases  in  the 
upper  atmosphere  ( e.g airglow)  (Munz  and  McFarland, 
1977).  The  final  irradiance  spectrum  is  strongly  red-shifted. 
Whereas  the  spectral  shift  from  moonlight  to  starlight  is 
minor  at  mesopelagic  depths,  it  is  quite  obvious  in  near- 
surface  waters  (Fig.  7A,  B),  particularly  in  blue,  oceanic 
waters.  Since  very  few  marine  species  are  known  to  have 
long-wavelength  sensitivity  at  scotopic  light  levels,  the  im¬ 
plications  of  the  spectral  shifts  at  these  wavelengths  are 
unknown.  However,  the  shifts  at  blue-green  wavelengths 
(450-500  nm)  are  also  substantial,  and  can  be  detected  by 
nearly  all  deep-sea  visual  systems.  Although  certain  coun- 
terilluminating  species  alter  the  spectra  of  their  emitted  light 
with  ambient  temperature  or  depth  (Young  and  Mencher, 
1980;  Young  and  Arnold,  1982;  Herring  et  al.,  1992), 
adaptations  to  the  spectral  shift  caused  by  the  presence  or 
absence  of  the  moon  are  unknown. 

Implications  for  counterillumination 

The  clarity  of  the  water  and  the  spectral  variation  due  to 
depth  and  the  presence  or  absence  of  the  moon  have  several 
important  implications  for  counterilluminators.  First,  since 
it  is  unlikely  that  light  scattering  by  the  water  will  combine 
the  light  from  the  individual  photophores  into  an  even  light 
field,  an  animal  with  few,  widely  spaced  light  organs  is  at  a 
disadvantage,  particularly  when  the  background  light  levels 
are  relatively  high.  Furthermore,  the  fewer  the  photophores, 
the  brighter  they  must  be  to  balance  out  the  unlit  regions  of 
the  ventral  surface.  In  this  study,  the  photophores  of  C. 
maderensis  had  to  be  175%  brighter  than  the  background 
radiance,  whereas  the  more  finely  distributed  photophores 
of  A.  veranyi  had  to  be  only  34%  brighter.  For  this  reason, 
a  counterilluminator  viewed  by  a  high-resolution  eye  will 
appear  as  a  signal  both  brighter  and  darker  than  the  back¬ 
ground  ( i.e .,  bright  photophores  on  a  silhouetted  body).  This 
may  explain  why  shallower  species  generally  have  more 
finely  spaced  photophores,  since  acute  vision  is  only  possi¬ 
ble  at  higher  levels  of  illumination  (Widder,  1999;  Warrant 
and  Locket,  2004). 

A  second  important  implication  of  this  study  is  that 
counterillumination  is  potentially  more  successful  at  shal¬ 
lower  depths.  Due  to  the  greater  contrast  attenuation  at 
shallow  depths,  any  mismatch  with  the  background  is  de¬ 
tectable  at  only  5%-20%  of  the  distance  at  which  the  same 
mismatch  would  be  detectable  in  deeper  waters.  This  in¬ 
crease  in  contrast  attenuation  may  offset  the  disadvantage 
due  to  the  variable  spectra  and  angular  distribution  found 
near  the  surface. 

Finally,  because  contrast  attenuation  is  relatively  constant 
over  a  wide  range  of  wavelengths  (Fig.  7D),  and  because 
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contrast  sensitivity  decreases  rather  slowly  with  decreasing 
illumination  (Warrant,  1999),  a  counterilluminator  ideally 
must  match  the  downwelling  spectrum  from  about  450  to 
520  nm  at  depth  and  over  a  somewhat  greater  wavelength 
range  near  the  surface.  However,  a  survey  of  published 
photophore  spectra  shows  that  this  is  not  the  case  (Fig.  8). 
Photophores  involved  in  counterillumination  do  have  spec¬ 
tral  characteristics  different  from  those  used  for  other  pur¬ 
poses,  but  the  pattern  is  far  from  intuitive.  In  fish,  counter- 
illuminating  photophores  produce  light  with  roughly  the 
same  peak  wavelength  (but  with  a  narrower  spectrum)  as 


those  of  non-counterilluminating  photophores.  In  decapods, 
the  peak  is  red-shifted  and  the  spectrum  narrower  in  coun- 
terilluminating  versus  non-counterilluminating  photo¬ 
phores.  In  squid  and  a  few  fish,  counterilluminating  photo¬ 
phores  emit  light  at  a  longer  (and  occasionally  shorter) 
wavelength  than  the  non-counterilluminating  photophores. 
Interestingly,  the  spectra  of  the  counterilluminators,  despite 
being  quite  clustered  (suggesting  natural  selection),  seldom 
match  the  downwelling  spectrum.  Some  are  10-20  nm  too 
blue,  and  others  are  10-30  nm  too  red.  This  suggests  that 
they  may  be  visible  to  predators  whose  color  discrimination 
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Figure  8.  Peak  wavelength  vs.  emission  width  (full-width  half-max,  FWHM)  for  the  light  emissions  from 
various  species.  White  symbols  denote  photophores  involved  in  counterillumination.  Black  symbols  denote 
photophores  and  other  luminous  sources  used  for  other  tasks  (warning,  luring,  etc  ).  The  outliers  among  the 
counterilluminators  are  xAbraliopsis  falco  (enoploteuthid  squid),  2Teuthowenia  megalops  (cranchiid  squid), 
*Isistius  hrasiliensis  (cookie- cutter  shark).  White  line  shows  peak  wavelength  and  FWHM  for  downwelling  light 
as  a  function  of  depth  (depth  intervals  are  10  m  down  to  100  m,  and  then  100  m  for  depths  down  to  800  m).  Bar 
chart  is  a  histogram  of  visual  pigment  absorption  peaks  for  deep-sea  fish  eyes  known  to  have  multiple  pigments 
in  the  blue-green  portion  of  the  spectrum  (data  from  Douglas  et  al.,  1998).  The  white  bars  are  the  short- 
wavelength  pigments;  the  black  bars  are  the  long-wavelength  pigments.  The  gray  symbols  show  the  cut-off 
wavelengths  for  the  filters  in  the  lenses  of  seven  species  of  deep-sea  fish  (data  from  Douglas  and  Thorpe,  1992). 
Photophore  spectral  data  taken  from  Nicol  (1960),  Swift  et  al.  (1973,  1977),  Biggley  et  al.  (1981),  Herring 
(1983),  Denton  et  al  (1985),  Widder  et  al.  (1983),  Herring  et  al.  (1992,  1993),  and  Haddock  and  Case  (1999). 
Because  the  published  bioluminescence  spectra  are  calibrated  in  energy  units  (instead  of  quanta),  the  down¬ 
welling  light  curve  is  calibrated  in  relative  energy  units. 
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at  blue  and  green  wavelengths  is  good  owing  to  multiple 
visual  pigments  or  ocular  filters. 

Although  about  15  deep-sea  species  (including  fish,  alci- 
opid  polychaetes,  oplophorid  shrimp,  and  the  squid  genus 
Watasenia)  are  known  to  have  multiple  visual  pigments, 
most  deep-sea  species  apparently  are  monochromatic,  with 
a  relatively  flat  spectral  response  near  the  maximum  wave¬ 
length  (due  to  the  extreme  thickness  of  their  photoreceptors) 
(Douglas  et  al,  1998).  Therefore  the  spectral  mismatches 
seen  in  Figure  8  may  not  be  detectable  as  a  color  shift  by 
many  predators  (excepting  those  possessing  color  filters). 
For  a  counterilluminator  facing  a  color-blind  predator,  the 
relevant  issue  is  that  the  light  emitted  from  the  photophores 
is  attenuated  more  quickly  than  the  downwelling  light,  due 
to  higher  absorption  at  non- peak  wavelengths.  Therefore, 
even  if  the  emitted  light  perfectly  matches  the  background 
intensity  at  one  distance,  the  counterillumination  will  be¬ 
come  darker  than  the  background  at  a  greater  distance.  The 
difference  between  the  attenuation  coefficients  at  470  nm 
and  480  nm  is  small.  Therefore,  this  issue  is  likely  to  be 
insignificant  for  the  krill,  fish,  and  decapods  whose  photo¬ 
phores  emit  at  470  nm.  The  close  clustering  of  the  spectra  of 
these  species  remains  puzzling,  but  may  be  an  evolutionary 
strategy  to  prevent  predators  from  developing  a  species- 
specific  search  image.  This  is  analogous  to  the  “selfish  herd” 
effect,  in  which  identical  individuals  in  large  aggregations 
reduce  their  chance  of  predation  (Hamilton,  1971;  Bond  and 
Al  Kamil,  2002). 

The  light  emitted  by  squid  photophores  that  peaks  at  510 
nm  will  be  attenuated  significantly  more  quickly  than  the 
downwelling  light,  potentially  leading  to  the  detection  of 
the  squid,  but  these  measured  spectra  may  not  be  represen¬ 
tative  of  the  natural  spectra.  As  mentioned  above,  certain 
squid  can  change  the  spectrum  of  their  counterillumination 
depending  on  temperature.  Since  the  spectral  measurements 
were  not  done  in  situ  and  often  required  a  fair  bit  of 
manipulation,  the  squid  may  have  been  at  a  higher  temper¬ 
ature  and  thus  produced  light  to  match  shallower  and  there¬ 
fore  greener  water.  In  fact,  the  published  spectra  of  all 
counterilluminators  must  be  treated  with  some  caution  be¬ 
cause  very  few  of  the  animals  were  measured  while  they 
were  passively  counterilluminating,  but  instead  were  being 
manually  stimulated  to  emit  light.  Because  manual  stimu¬ 
lation  tends  to  turn  on  all  available  photophores  in  an 
attention-getting  signal  that  is  assumed  to  act  as  a  “burglar 
alarm”  (Widder,  1999),  the  measured  spectrum  may  include 
light  from  photophores  that  are  not  active  during  counter¬ 
illumination,  altering  the  spectrum. 

Effects  of  visual  resolution  and  color  discrimination 
on  perceived  image 

Although  the  range  of  water  types  commonly  inhabited 
by  counterilluminators  had  little  effect  on  their  visibility, 


the  range  of  visual  acuities  of  potential  predators  had  a 
dramatic  effect.  Because  light  scattering  by  the  water  had 
little  effect  on  the  appearance  of  the  counterillummation 
signal,  acute  vision  can  detect  the  individual  photophores 
and  is  therefore  highly  advantageous.  Indeed,  many  deep- 
sea  species  are  known  to  have  far  greater  resolution  (~10X) 
in  the  dorsal  viewing  region  than  in  other  directions  (Collin 
et  al.,  1997;  Wagner  et  al.,  1998;  Land,  1999;  Warrant  and 
Locket,  2004).  For  example,  although  the  spatial  resolution 
for  upward  viewing  in  the  hatchet  fish  Argyropelecus  ac - 
uleatus  is  0.11°  (see  Materials  and  Methods),  the  spatial 
resolution  over  the  rest  of  visual  field  is  0.7°-1.7°  (Collin  et 
al.,  1997).  In  contrast,  the  myctophid  Lampanyctus  festivus, 
which  has  lateral-viewing  eyes,  has  a  relatively  constant  and 
low  visual  acuity  (0.5°)  over  the  entire  visual  field  (Wagner 
et  al.,  1998).  Because  this  increased  spatial  resolution  de¬ 
creases  sensitivity  (and  hence  ability  to  detect  contrast),  it 
has  an  associated  cost.  Warrant  and  Locket  (2004)  analyzed 
the  benefits  and  costs  of  high  spatial  resolution  as  a  function 
of  what  is  being  viewed;  they  determined  that  high  spatial 
resolution  should  be  selected  for  in  eyes  that  search  over¬ 
head  for  small,  silhouetted  objects.  While  they  do  not  ex¬ 
plicitly  consider  the  spatial  pattern  of  counterillumination, 
the  same  principles  apply. 

The  high  spectral  variation  of  the  background  light  and 
the  spectral  mismatches  seen  in  Figure  8  suggest  that  good 
color  discrimination  in  the  blue-green  would  be  extremely 
advantageous.  As  mentioned  above,  certain  deep-sea  spe¬ 
cies  probably  have  good  color  discrimination  at  blue  and 
green  wavelengths.  Indeed,  the  peaks  of  these  pigments 
seem  to  support  the  hypothesis  of  Douglas  et  al.  (1998)  that 
one  pigment  matches  the  counterilluminator’s  spectrum  and 
one  matches  the  downwelling  light  (Fig.  8).  In  addition, 
certain  species  with  only  one  visual  pigment  have  multi- 
banked  retinas.  The  filtering  of  the  light  by  the  vitread  banks 
alters  the  spectrum  of  the  light  reaching  the  sclerad  banks, 
changing  their  sensitivity  and  theoretically  allowing  for 
color  discrimination  (Denton  and  Locket,  1989).  Finally, 
the  lenses  of  certain  deep-sea  species  have  yellow  filters  that 
can  also  increase  the  contrast  of  a  counterilluminator  against 
the  background  (Munz,  1976;  Douglas  and  Thorpe,  1992). 

Conclusions 

Although  counterillumination  is  a  ubiquitous  and  suc¬ 
cessful  cryptic  strategy,  the  clarity  of  the  water  implies  that 
the  camouflage  can  be  broken  by  species  with  acute  vision 
at  longer  distances  than  anticipated.  In  addition,  the  back¬ 
ground  to  be  matched  depends  not  only  on  depth,  but  also 
on  the  source  of  nocturnal  illumination.  While  spectral 
variation  is  greatest  near  the  surface,  contrast  attenuation  is 
also  greatest.  These  results  suggest  several  fruitful  avenues 
for  future  research,  including  further  analysis  of  the  con¬ 
flicting  constraints  of  visual  sensitivity  and  spatial  resolu- 
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tion,  a  determination  of  how  counterilluminators  that  can 
change  spectral  emissions  choose  the  correct  spectrum  (de¬ 
spite  being  color-blind),  and  investigation  of  a  possible 
relationship  between  the  resolution  of  ventral  photophore 
patterns  and  the  acuity  of  potential  predators. 
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Abstract.  Artificial  radiance  sets  were  used  as  inputs  to  Multi-layer  Perceptron 
and  multilinear  regression  algorithms  to  study  their  retrieval  capabilities  for 
optically  active  constituents  in  sea  water.  The  radiative  transfer  model 
Hydrolight  was  used  to  produce  18,000  artificial  reflectance  spectra  representing 
various  case  1  and  case  2  water  conditions.  The  remote  sensing  reflectances  were 
generated  at  the  Medium  Resolution  Imaging  Spectrometer  (MERIS)  wave¬ 
bands  412,  442,  490,  510,  560,  620,  665  and  682  nm  from  randomly  generated 
triplet  combinations  of  chlorophyll  a,  non-chlorophyllous  particles  and  CDOM 
(Coloured  Dissolved  Organic  Matter)  concentrations.  These  reflectances  were 
contaminated  with  different  noise  terms,  before  they  were  used  to  assess  the 
performance  of  multilayer  perceptron  and  multilinear  regression  algorithms.  The 
potential  of  both  algorithms  for  retrieving  optically  active  constituents  was 
demonstrated  with  the  neural  network  showing  more  accurate  results  for  case  2 
scenarios. 


1.  Introduction 

With  the  deployment  of  the  Medium  Resolution  Imaging  Spectrometer 
(MERIS)  instrument  on  board  the  Envisat  satellite  from  the  European  Space 
Agency,  a  need  has  arisen  for  more  complex  ocean  colour  algorithms  that  cater  for 
the  larger  selection  of  bands  that  MERIS  has  to  offer  compared  to  similar 
instruments.  This  should  result  in  more  accurate  results  for  coastal  regions  than 
from  previous  ocean  colour  instruments.  Coastal  regions  are  very  dynamic 
environments  where  conditions  vary  over  short  time  and  spatial  scales.  Tidal 
currents  and  river-runoff  carrying  coloured  dissolved  organic  matter  (CDOM)  and 
suspended  sediment  into  the  oceans  influence  ocean  colour  in  these  waters  (Gould 
and  Amone  1997).  Coastal  waters  are  typically  classified  as  case  2,  signifying  that 
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chlorophyll  a  is  not  the  only  significant  influence  on  the  colour  and  there  are  other 
optically  active  substances  not  correlated  with  the  chlorophyll  a  (Mobley  1994).  In 
such  waters,  the  concentration  of  the  optically  active  constituents  needs  to  be 
estimated  from  the  surface  colour  spectra.  Numerous  algorithms  have  been 
developed  to  undertake  this  process  and  recently  some  based  on  multi-layer 
perceptron  (MLP)  neural  networks  have  shown  promising  capabilities  (Gross  et  al 
2000).  The  purpose  of  this  research  letter  is  to  compare  the  ocean  colour  parameter 
estimation  performance  of  neural  networks  and  multilinear  regression  algorithms 
presented  with  ocean  colour  spectra  obtained  from  the  ocean  colour  model 
Hydrolight. 

2.  Artificial  radiance  spectra  generation 

The  data  used  for  this  research  were  artificially  generated  as  an  alternative  to 
the  use  of  in  situ  measurements  of  water  parameters  matched  to  a  concurrent 
satellite  overpass.  This  approach  was  adopted  in  order  to  have  control  over  the 
amount  of  noise  in  the  spectral  sets.  In  total  18,000  remote  sensing  reflectance 
spectra  were  produced;  9000  for  case  1  and  9000  for  case  2  conditions.  This  was 
carried  out  in  three  steps  described  in  the  following  three  sections.  Section  2.4 
describes  the  noise  terms  added  to  the  data. 

2.1.  Generation  of  optically  active  parameter  concentrations 

Chlorophyll  a,  non-chlorophyllous  matter  and  CDOM  can  be  expressed  as  the 
optically  active  parameters  (OAPs)  in  sea  water.  Their  parameterisations  are  the 
chlorophyll  a  mass  to  volume  concentration  C  (mg  m“3),  the  scattering  coefficient 
of  non-chlorophyllous  particles  at  550  nm  X  (m-1)  and  the  absorption  coefficient  of 
yellow  substance  at  440  nm  Y  (m-1).  A  combination  of  C,  X  and  Y  is  used  to 
produce  one  spectrum.  A  random  number  generator  has  been  used  to  create  the 
OAP  combinations  according  to  the  statistical  distribution  parameters  shown  in 
table  1,  where  p  and  a  represent  the  distribution  mean  and  standard  deviation 
respectively.  The  values  for  case  1  water  in  table  1  have  been  based  on  the  analysis 
of  a  CZCS  (coastal  zone  colour  scanner)  global  composite  pigment  map  histogram 
and  the  case  2  water  values  reflect  typical  case  2  concentration  ranges  (Cipollini 
1996). 

2.2.  The  bio-optical  model 

To  evaluate  the  scattering  and  absorption  values  in  the  water  related  to  each 
combination  of  OAP  values,  a  biooptical  model  was  used.  Total  absorption  a  was 
partitioned  into  the  separate  contributions  of  pure  seawater  aWi  phytoplankton  aph. 


Table  1.  Statistical  properties  of  OAP  generation. 


Case  I  waters 

Case  II  waters 

PLog 

&  Log 

P-Log 

°Log 

Log(C) 

-0.86 

0.3 

0.0 

0.5 

Log(X) 

-1.21 

0.3 

0.0 

0.5 

Log(Y) 

-1.75 

0.3 

-0.5 

0.5 

PLog(C).Log(X) 

PLog(C).Log(Y) 

PLog(C).Log(X) 

PLog(C).Log(Y) 

correlation 

0.8 

0.8 

0.5 

0.5 
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sediment  ased  and  CDOM  ay.  The  scattering  coefficient  was  partitioned  into  the 
separate  contributions  of  water  bW)  phytoplankton  bph  and  sediment  scattering  bsed. 

The  biooptical  model  was  used  to  obtain  the  total  absorption  and  scattering 
coefficients  for  the  eight  MERIS  bands.  Once  these  had  been  calculated,  they  were 
used  as  the  inputs  to  the  radiative  transfer  model  that  converted  the  absorption  and 
scattering  coefficients  into  the  remote  sensing  reflectances.  The  detailed  biooptical 
model  can  be  found  in  Cipollini  et  al.  (2001).  The  phytoplankton  scattering 
coefficient  included  a  fluorescence  peak  directly  proportional  to  the  C  concentration 
which,  in  this  letter,  has  been  amended  slightly  to  avoid  a  fluorescence  peak  solely 
dependent  on  the  chlorophyll  a  concentration.  Expression  (1)  shows  the  amended 
scattering  coefficient  for  phytoplankton  based  on  the  original  expression  given  by 
Cipollini  et  al  (2001). 


(i) 


The  amendment  consists  of  a  random  variable  R  controlling  the  height  of  the 
fluorescence  peak  that  was  originally  modelled  by  superimposing  a  Gaussian  peak 
onto  the  phytoplankton  scattering  spectrum  commonly  modelled  as  a  flat  line. 

2.3.  Radiative  transfer  model 

The  absorption  and  scattering  values  computed  by  the  biooptical  model  were 
used  as  an  input  to  the  radiative  transfer  model  Hydrolight  (Mobley,  1995)  that 
generated  the  remotely  sensed  reflectance  for  each  OAP  combination.  Two  datasets 
for  case  1  and  case  2  waters,  respectively,  were  created  each  totalling  9000  spectra 
and  their  corresponding  OAP  combinations. 

2.4.  Noise  terms  added  to  datasets 

Three  different  noise  terms  were  added  to  the  datasets  to  represent  noise  due  to 
the  remote  sensing  instrument,  an  in  situ  measurement  noise  inherent  to  any 
instrument  calibration  datasets  and  an  atmospheric  correction  error  residual  to 
ocean  colour  products  after  the  atmospheric  contributions  have  been  removed.  The 
instrument  noise  is  a  Gaussian  noise  within  ±2%  of  the  reflectance  value  in  each 
band  and  the  in  situ  noise  is  a  Gaussian  noise  within  ±40%  of  the  concentration 
value  of  each  OAP.  In  both  cases  the  percentages  represent  the  extremes  the 
noiselevels  can  reach.  The  atmospheric  correction  error  was  simulated  by  taking 
three  reflectance  values  at  the  bands  442,  490  and  560  nm,  suggested  by  Antoine 
and  Morel  (1999)  to  be  the  largest  allowable  errors  for  an  atmospheric  correction, 
and  using  these  values  as  negative  and  positive  limits  for  a  Gaussian  error 
distribution.  The  values  for  the  remaining  bands  were  found  by  linear  interpolation 
and  extrapolation  so  that  the  maximum  error  spectrum  shown  in  table  2  was 
obtained. 

3.  Algorithms 

The  MLP  neural  network  used  for  the  case  1  data  set  had  8  input  nodes,  4 
nodes  in  the  first  hidden  layer,  3  in  the  next  and  3  output  nodes  for  each  OAP 
derived.  A  sigmoid  nodal  transfer  function  was  used  for  both  case  1  and  case  2 
data.  The  case  2  MLP  was  chosen  to  be  a  bit  more  complex  so  that  it  could  deal 
with  the  greater  spread  of  the  case  2  data.  It  consisted  of  8  input  nodes,  6  nodes  in 
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Table  2.  Atmospheric  noise  boundaries. 


Wavelength  (nm) 

Reflectance  +/— 

412 

0.003 

442 

0.002 

490 

0.0005 

510 

0.0004 

560 

0.0002 

620 

0.0002 

665 

0.0002 

682 

0.0002 

the  first  hidden  layer,  three  nodes  in  the  second  hidden  layer  and  3  output  nodes. 
These  architectures  were  determined,  by  a  systematic  procedure  of  reducing 
complex  nets  to  simpler  ones  without  a  loss  of  performance.  In  that  procedure, 
nodes  were  removed  from  the  MLPs  until  their  performance  in  terms  of  relative 
RMS  errors  and  determination  coefficients,  started  deteriorating,  i.e.  the  errors 
increased  and  the  determination  coefficients  decreased.  A  similar  procedure  was 
carried  out  to  find  optimum  training  and  query  set  sizes  as  well  as  the  optimum 
number  of  iterations  before  over  fitting  becomes  an  issue.  For  both  case  1  and  case 
2  the  training  and  query  set  sizes  were  6000  and  3000  spectra,  respectively,  as  well 
as  a  training  time  of  6000  iterations. 

The  multilinear  regression  as  shown  in  (2)  was  carried  out  on  the  6000  spectra 
used  for  the  training  set,  and  the  3000  spectra  for  the  test  set  were  used  to  assess  the 
performance  of  the  regression. 

N 

log(C)  =  a0+  ^a.log(-R,)  (2) 

1=1 


4.  Results  for  case  1  and  case  2  waters 

The  capacity  of  the  algorithms’  outputs  to  match  the  original  spectra-generating 
OAPs  is  presented  in  terms  of  the  relative  root-mean-square  error 

-2 


/|  n 

relative  rms  =  4 

\  nk=  l 


tk-ok 


tk 


indicating  the  precision  of  the  predicted  values,  and  the  coefficient  of  determination 

R2: 


R2  =  l~J^’  Where  SSE  =  £  O'. -h?  and  SST  =  (£ yfj  - 

SSE  =  sum  of  squares  error,  SST  =  total  sum  of  squares 
y,  =  datapoint,  y,  =  adjunct  point  on  trendline 
t  =  target  value,  o  =  output  value  from  neural  network 


The  results  of  all  algorithms  for  the  case  1  water  set  can  be  seen  in  table  3a  &  b. 
For  case  1,  the  neural  network  and  multilinear  algorithm  perform  similarly  with 
large  R2  values  and  small  RMS  errors.  The  MLP  is  showing  an  overall  more 
accurate  performance  despite  a  slightly  larger  RMS  error  for  the  estimation  of  C. 
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Algorithm 

C  rRMS 

X  rRMS 

Y  rRMS 

C  R 2 

X  R2 

Y  R2 

MLP 

0.220 

0.208 

0.178 

0.902 

0.912 

0.903 

ML 

0.197 

0.219 

0.215 

0.869 

0.881 

0.855 

Table  3(b).  rel.  RMS 

errors  and  R 2 

for  case  2. 

Algorithm 

C  rRMS 

X  rRMS 

Y  rRMS 

C  R2 

X  R2 

Y  R2 

MLP 

0.539 

0.418 

0.226 

0.651 

0.953 

0.920 

ML 

0.909 

0.240 

0.290 

0.557 

0.938 

0.825 

For  case  2  the  MLP  is  most  accurate  for  the  estimation  of  C  and  Y  whereas  the  ML 
algorithm  estimates  X  more  precisely  indicated  by  the  smaller  RMS  error. 

5.  Conclusion 

Whereas  both  algorithms  show  similar  performance  levels  for  case  1  data,  the 
MLP  shows  a  more  significant  potential  for  the  retrieval  of  OAPs  from  reflectance 
spectra  of  the  case  2  dataset.  Hence  the  authors  propose  that  for  ocean  colour 
algorithm  development,  emphasis  should  be  placed  on  the  use  of  MLP  neural 
networks.  However,  the  relatively  inaccurate  results  also  evident  for  the  MLP, 
during  the  C  retrieval  of  case  2  water  parameters  indicate  the  difficulty  associated 
with  these  waters.  One  of  the  main  challenges  for  ocean  colour  algorithm 
development  therefore  is  the  creation  of  algorithms  capable  of  generating  accurate 
ocean  colour  products  from  case  2  regional  imagery.  The  case  2  data  generated  for 
this  letter  are  generic  representing  a  wide  variety  of  conditions.  The  great  variability 
that  exists  between  different  case  2  regions  in  contrast  suggests  that  algorithms 
should  be  developed  using  regionally  specific  models. 
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[l]  It  is  shown  that  the  in- water,  shape  factor  formulation  of  the  radiative  transfer 
equation  (RTE)  (1)  yields  exact  in-air  expressions  for  the  remote  sensing  reflectance  RTS 
and  the  equivalent  remotely  sensed  reflectance  RSRa  and  (2)  can  be  configured  for 
inherent  optical  property  (IOP)  retrievals  using  standard  linear  matrix  inversion  methods. 

Inversion  of  the  shape  factor  RTE  is  exact  in  the  sense  that  no  approximations  are  made  to 
the  RTE.  Thus  errors  in  retrieved  IOPs  are  produced  only  by  uncertainties  in  (1)  the 
models  for  the  shape  factors  and  related  quantities  and  (2)  the  IOP  models  required  for 
inversion.  Hydrolight  radiative  transfer  calculations  are  used  to  derive  analytical  models 
for  the  necessary  backscattering  shape  factor,  radiance  shape  factor,  fractional  forward 
scattering  coefficient,  ratio  of  air-to-water  mean  cosines,  and  diffuse  attenuation 
coefficient  for  in-water  upwelling  radiance.  These  models  predict  the  various  shape  factors 
with  accuracies  ranging  typically  from  2  to  20%.  Using  the  modeled  shape  factors  the 
in-air  remotely  sensed  reflectance  RSRa  can  be  predicted  to  within  20%  of  the  correct 
(Hydrolight-computed)  values  96%  of  the  time  (or  ±0.0005  sr"1  86%  of  the  time)  for  the 
synthetic  data  used  to  determine  the  shape  factor  models.  Inversion  of  this  shape 
factor  RTE  using  field  data  is  a  comprehensive  study  to  be  published  in  a  later 
paper.  INDEX  TERMS:  4552  Oceanography:  Physical:  Ocean  optics;  4847  Oceanography:  Biological  and 
Chemical:  Optics;  4275  Oceanography:  General:  Remote  sensing  and  electromagnetic  processes  (0689);  4842 
Oceanography:  Biological  and  Chemical:  Modeling;  KEYWORDS:  remote  sensing,  optical  oceanography, 
inverse  modeling,  radiative  transfer  theory 
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1.  Introduction 

[2]  Semianalytic  radiance  models  [Gordon  et  al. ,  1988; 
Morel  and  Gentili,  1996]  can  be  readily  inverted  by  linear 
matrix  methods  [Hoge  etal .,  1999a,  1999b,  2001]  to  provide 
oceanic  inherent  optical  properties  (IOPs).  Such  inversions 
are  well  conditioned  [Hoge  and  Lyon ,  1996]  and  promise  a 
powerful  method  of  simultaneously  retrieving  constituent 
absorption  and  backscattering  coefficients  in  the  upper 
surface  layer  of  the  world’s  oceans  using  satellite  data  [Hoge 
et  al ,  2001;  Hoge  and  Lyon ,  2002].  However,  semianalytic 
radiance  models  (1)  do  not  provide  an  exact  framework 
to  account  for  all  possible  environmental  and  viewing 
conditions  [Weidemann  et  al.,  1995]  and  (2)  contain  fixed 
constants  that  both  obscure  insight  into  the  physical  radiative 
transfer  processes  and  limit  their  flexibility. 
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[3]  The  radiative  transfer  equation  (RTE)  can  provide 
exact  inverse  solutions,  but  the  RTE  is  not  easily  inverted 
for  many  remote  sensing  situations  [Zaneveld,  1995]. 
Therefore  a  specific  form  of  the  RTE  inversion  is  investi¬ 
gated,  namely  a  modified  version  of  the  shape  factor 
formulation  of  Zaneveld  [1995].  Some  of  the  motivation 
for  the  work  herein  comes  from  the  distinct  need  for  highly 
accurate  methods  to  retrieve  the  absorption  coefficients  of 
the  chlorophyll  accessory  pigment  phycoerythrin  [Hoge  et 
al.,  1999b].  To  this  end  the  absorption  coefficients  of 
chlorophyll  and  chromophoric  dissolved  organic  matter 
(CDOM)  must  be  accurately  retrieved;  otherwise,  weaker 
absorbing  constituents  (such  as  phycoerythrin)  will  be 
obscured. 

[4]  In  this  paper  (1)  the  shape  factor  form  of  the  RTE  is 
shown  to  be  readily  configured  into  linear  form  for  simul¬ 
taneous  retrieval  of  oceanic  IOPs  using  standard  matrix 
methods;  (2)  the  RTE  inversion  is  derived  for  the  principal 
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“big  three”  IOPs,  namely  the  phytoplankton  absorption 
coefficient,  the  CDOM  +  detritus  absorption  coefficient, 
and  the  total  constituent  backscattering  coefficient;  (3)  shape 
factor  and  related  models  required  for  the  inversion  are 
developed  for  backscattering  and  radiance  shape  factors,  the 
diffuse  attenuation  coefficient  for  upwelling  radiance,  the 
ratio  of  average  cosines  of  the  air  and  water  downwelling 
irradiances,  and  the  fractional  forward  scattering  coefficient; 
and  (4)  propagation  of  errors  into  the  IOP  state  vector 
resulting  from  errors  in  the  data-model  matrix  and  hydro- 
spheric  vector  as  well  as  shape  factor  and  related  models  are 
assessed. 

[5]  Our  ultimate  objective  is  to  determine  if  the  shape 
factor  RTE  matrix  inversion  methodology  will  result  in 
accurate  algorithms  for  application  to  satellite  ocean  color 
data  This  paper  presents  the  underlying  shape  factor  RTE 
theory  and  develops  the  needed  models  for  the  shape  factors 
and  related  quantities,  while  future  work  will  describe 
comprehensive  studies  of  the  shape  factor  RTE  inversion 
of  synthetic  and  real  data. 


2.  Shape  Factor  Form  of  the  Radiative 
Transfer  Equation 

[6]  Establish  a  Cartesian  coordinate  system  with  +z  axis 
vertically  downward  into  the  ocean  and  x  and  y  axes  lying 
within  the  atmosphere-ocean  boundary.  In  a  plane  parallel 
medium  without  internal  sources  or  inelastic  scattenng,  the 
radiative  transfer  equation  is 

cosO^^—1 —  =  —  c(z)Z,(0,  <j>,z) 
dz 

2k  k 

+  J  yp(e,(t);0',<|>';*)L(0',<t),,r)sin0'  d&  dtf. 
0  0 

(1) 

(See  notation  section  for  definition  of  symbols.)  Zaneveld 
[1995,  1982]  showed  that  equation  (1)  can  be  rewritten  in 
terms  of  the  in-water  remotely  sensed  reflectance  (RSR)  as 


RSR  = 


/»(0.  <!>.*) 


Mf) 

2t\ 


-cose*(e,<t >,z)  +c(z)  -/i(e,<MMz)  ’ 

(2) 


where  the  dimensionless  backscattering  shape  factor  fb(Q, 
(j),  z)  is  given  by 


2k  k/2 
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0  0 
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the  dimensionless  radiance  shape  factor  ^(0,  4>,  z)  is  given 
by 
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and  the  diffuse  attenuation  coefficient  for  upwelling 
radiance  &(0,  <)>,  z),  with  units  of  m-1,  is  given  by 


1 


</Lu(0,  4>,  zr) 


£u(6,<l>,  z)  dz 


(5) 


Equation  (2)  is  Zaneveld' s  [1995]  equation  (7)  and  is  exact 
because  it  is  simply  a  restatement  of  the  RTE  (1)  for  upward 
directions  using  definitions  (3)-(5).  Subscripts  d  and  u 
appended  to  the  radiance  L  explicitly  remind  us  that  the 
radiance  in  equation  (3)  is  downwelling  (0  <  0  <  tt/2), 
whereas  the  radiance  in  equations  (4)  and  (5)  is  upwelling 
(tt/2  <  0  <  tt)-  (The  iso-  superscript  is  discussed  below.) 

[7]  The  numerator  of  the  fb  shape  factor  in  equation  (3) 
shows  how  much  downwelling  radiance  is  scattered  upward 
into  direction  (0,  <j>).  The  denominator  is  the  same  quantity 
evaluated  for  the  special  case  of  an  isotropic  volume 
scattering  function  (in  which  case  (3  =  2b^/di().  Thus  the 
fb  shape  factor  is  a  measure  of  how  much  the  actual  phase 
function  differs  from  a  constant  over  the  backscattering 
directions.  Similarly,  the  numerator  of  fL  in  equation  (4) 
shows  how  much  the  upwelling  radiance  is  forward  scat¬ 
tered  into  direction  (0,  <j>).  The  denominator  is  the  same 
quantity  evaluated  for  the  special  case  of  an  isotropic 
upwelling  radiance  distribution  whose  magnitude  is  E™ 
and  for  the  special  case  of  an  isotropic  volume  scattering 
function.  Clearly,  these  shape  factors  depend  both  on  the 
IOPs  (namely  on  the  volume-scattering  function,  in  this 
case)  and  on  the  ambient  radiance  distribution,  as  does  the 
diffuse  attenuation  coefficient  of  equation  (5).  These  quan¬ 
tities  therefore  are  unknown  terms  in  equation  (2)  if  equa¬ 
tion  (2)  is  to  be  inverted  to  obtain  the  IOPs  a  and  bh  from 
measured  upwelling  radiances  and  downwelling  irradiances. 
The  fact  that  shape  factors  are  unknown  prevents  the  RTE  in 
equation  (2)  from  being  inverted  unless  further  assumptions 
are  made  about  the  values  of  the  shape  factors.  Modeling 
these  unknowns  in  terms  of  known  quantities  is  the  major 
focus  of  this  paper. 

[8]  Equations  (1 ) — (5)  are  valid  at  any  depth  within  an 
arbitrarily  stratified  water  column,  but  the  specific  interest 
herein  is  remote  sensing  of  near-surface  water  IOPs  There¬ 
fore  one  needs  to  relate  the  quantities  in  equations  (2)- (5), 
when  evaluated  just  beneath  the  mean  sea  surface,  to 
quantities  in  air  just  above  the  sea  surface,  which  can  be 
deduced  via  in-air  remote  sensing  techniques.  Equation  (2) 
can  be  converted  into  a  form  suitable  for  above-water 
remote  sensing  applications  as  follows.  The  ^-squared  law 
for  radiance  transmittance  across  a  boundary  between  two 
media  [Mobley,  1994,  equation  (4.21)]  can  be  used  to 
convert  the  in-water  upwelling  radiance  just  beneath  the 
sea  surface,  Z,„(0,  (j>,  z  =  0),  to  the  water-leaving  radiance  in¬ 
air  just  above  the  sea  surface,  Lua(Qa,  4>): 

L„(0,<M  =  O)=^M0ai<t,).  (6) 

Subscript  a  denotes  values  in  air,  just  above  the  mean  sea 
surface;  depth  z  =  0  denotes  values  in  water,  just  beneath 
the  sea  surface.  The  in-air  polar  angle  0a  associated  with 
Euafiay  4>)  is  the  refracted  viewing  angle  above  the  sea 
surface  obtained  by  applying  Snell’s  law  to  the  in-water 
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angle  0.  The  downwelling  scalar  lrradiance  £od(z)  is 
converted  to  the  downwelling  plane  irradiance  £d(z)  via 
the  mean  cosine  of  the  downwelling  radiance,  Eod(z )  = 

£d(z)/ \x,d.  The  plane  irradiance  just  beneath  the  sea  surface 
can  be  related  to  the  in-air  value  via  [Mobley,  1994, 
equation  (7. 19)]  2sd(z  =  0)  =  EdJ/ (1  +  rR).  Combining  these 
results  gives  [Mobley,  1994,  equation  (10.27)] 


Call  Dl  the  radiance  derivative  term  because  it  is  a  measure 
of  the  depth  rate  of  change  of  the  upwelling  radiance,  as 
seen  in  equation  (5).  Define  the  second  term  in  the 
denominator  as 

MM*0)  =  M°)[i  -A(e,<M)].  (n) 


Z,„(e,(t>,  z  =  0)  _  nj(  1  -rR)\jd  L,,o(0fl,4>)  ,  . 

Eod{z  =  0)  1 1  Eda  U 

In  equation  (7),  define  M  =  [(t  7)/[(l  -  rR)rtJ]  For  a  wide 
range  of  sky  and  sea  surface  conditions  and  for  viewing 
directions  relevant  to  remote  sensing,  M  lies  in  the  range  of 
0.53  to  0.55  [Mobley,  1994;  Hoge  and  Lyon,  1996;  Hoge  et 
al,  1999a,  1999b;  Morel  and  Gentili,  1996].  Thus  Mean  be 
approximated  a sM«  0.54,  with  an  error  of  less  than  2%. 
Using  equation  (7)  and  partitioning  the  beam  attenuation 
coefficient  as  c(z)  =  a(z)  +  6/(z)  +  bb(z),  equation  (2)  becomes 


The  shape  factor^  varies  from  0.963  to  1.152  for  nadir 
viewing  [Weidemann  et  al,  1995;  Zaneveld,  1995];  thus  Bf 
ranges  from  0.0372yto  —0.152 bf,  which  is  a  small  fraction 
of  the  forward  scattering  coefficient  bf.  Therefore  call  Bf  the 
fractional  forward  scattering  coefficient.  Thus  bf  and  fL  are 
found  in  a  combination  in  which  one  (f£)  serves  to  reduce 
the  size  of  the  other  (bf).  Finally,  define  the  mean  cosine 
ratio  as 


(12) 


n  LUo(*aA) 

Kr*  ~  — ~Ed - 

M/fc(M,o)Mo)/[2*iv(o)] _ 

-*(0, 4),  0)  cos  0  +  6/(0) [1  -/L(0, 4),  0)]  +  a( 0)  -P  bb( 0)  ■ 

(8) 

Except  for  the  small  error  associated  with  the  assumed  value 
for  M,  equation  (8)  remains  an  exact  RTE  expression  for  the 
in-air  remote  sensing  reflectance,  Rrs ,  just  above  the  sea 
surface.  The  remote  sensing  reflectance  is  the  quantity  used 
as  the  basis  for  ocean  color  remote  sensing  by  the  Sea¬ 
viewing  Wide  Field-of-view  Sensor  (SeaWiFS)  [O’Reilly  et 
al.,  1998;  Hoge  et  al.,  2001;  Hoge  and  Lyon,  2002]  and 
airborne  systems  [Davis  et  al,  2002;  Hoge  et  al,  1999a, 
1999b]  systems.  Rxs  can  be  obtained  from  at-sensor  radiances 
after  atmospheric  correction;  for  our  purposes  here  it  is 
therefore  considered  known.  As  noted  by  Zaneveld  [  1995],  it 
is  desirable  to  use  RSRa,  the  in-air  value  of  RSR,  rather  than 
Rrs  because  the  scalar  irradiance  Eod  is  less  sensitive  to  solar 
zenith  angle  effects  than  is  the  plane  irradiance  £d.  Thus  use 
RSRa  =  \\jaRr*  to  rewrite  equation  (8)  as 


RSRa  = 


Eoda 


-*(0, 4>,  0)  cos 0  +  6/(0)[l  -/i(0, <j>, 0)]  +  a(0)  +  bb( 0)  ' 

(9) 


Using  definitions  (1 0)— (1 2),  equation  (9)  becomes 


RSR^,^)  = 


E0da 


Dl(6, <t>,  0)  +  5/(0,  <!>,  0)  +  a( 0)  4-  M0) ' 


(13) 


Equations  (8),  (9),  and  (13)  are  each  called  the  shape  factor 
form  of  the  RTE.  Equation  (13)  is  addressed  hereafter.  The 
ultimate  goal  is  to  use  equation  (13)  to  relate  the  unknown 
absorption  and  backscattering  coefficients  just  beneath  the 
sea  surface  to  the  known  remotely  sensed  reflectance  and 
other  known  quantities.  As  noted  above,  M  =  0.54. 
However,  the  four  quantities  fb,  R^,  DL,  and  Bf  (or, 
equivalently,  fb,  R ^  k,  and  6/(1  -  fL)  as  seen  in 
equation  (9)  and  for  brevity  call  all  of  these  quantities 
shape  factors)  are  unknown.  The  shape  factors  depend  in 
complicated  ways  on  the  water  column  IOPs,  environmental 
conditions  (sky  radiance  and  sea  state),  and  viewing 
geometry  (Sun  zenith  angle  and  viewing  direction).  In 
section  4  the  shape  factors  are  modeled,  so  that  they  too  can 
be  considered  known  in  equation  (13). 


3.  Linear  Form  of  the  Radiative  Transfer 
Equation  and  Its  Inversion 

[10]  The  m-air  RSRa  of  equation  (13)  immediately  yields 
the  fundamental  linear  form  of  the  RTE, 


The  simplicity  of  Zaneveld' s  [1995]  original  in-water  RSR 
formulation  remains  in  this  equation  for  RSRa,  except  for  M 
and  the  fw'prf  ratio  for  the  downwelling  light  field.  As  a 
practical  matter,  equation  (8)  is  presently  more  easily  applied 
to  oceanic  field  data  because  the  in-air  downwelling  plane 
irradiance  is  more  generally  available,  but  there  are  no 
instrumental  barriers  to  using  scalar  irradiance  as  in 
equation  (9). 

[9]  To  further  simplify  equation  (9)  for  later  use,  define 
the  first  term  in  the  denominator  as 

£>i(e,4>>0)  =  -*(0>4>,O)cos0.  (10) 


a(0)  +  6*(0)F  +  DL+R/  =  0, 


where 


F(0,4>,O)  =  1- 


M  — Jit - 

Luai^a)  4>) 

Soda 


(14) 


(15) 


V  is  called  the  backscattering  enhancement  factor. 

[11]  Next,  partition  the  total  absorption  coefficient  into 
contributions  by  pure  water,  phytoplankton,  and  CDOM 
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plus  detritus.  Similarly,  the  backscattering  coefficient  is 
written  as  the  sum  of  contributions  by  pure  seawater  and 
by  particulate  matter.  It  is  easy  to  show  [Hoge  and  Lyon , 
1996;  Hoge  et  al. ,  1999a,  1999b]  that  the  equation  describ¬ 
ing  the  desired  phytoplankton  absorption  coefficient  aph , 
CDOM  +  detritus  absorption  coefficient  adi  and  total 
constituent  backscattering  coefficient  bbt  resulting  from 
equation  (14)  is 

aPh{\)  4-  ad(\i)  +  bbt(\i)V  =  -aw(X,)  -  bbw(\)V  -  Di  —  Bf. 

(16) 

The  wavelength  dependency  of  the  IOPs  is  now  shown 
explicitly,  while  the  depth  and  angular  dependencies  have 
been  suppressed  for  clarity.  Note  that  the  observed  water- 
leaving  radiances  Lua  occur  on  both  sides  of  the  equation 
(within  V).  The  pure  water  absorption  a„  is  known  from 
Pope  and  Fry  [1997],  and  the  water  backscattering 
coefficient  bbw  is  given  by  Smith  and  Baker  [1981].  The 
right-hand  side  of  equation  (16)  is  therefore  known,  given 
the  shape  factors  and  a  measurement  of  RSRa(\).  This 
linear  form  of  the  RTE  is  still  exact  in  the  sense  that  no 
approximations  have  been  made  to  the  RTE,  but  clearly,  the 
IOP  retrieval  accuracy  will  be  determined  by  the  accuracy 
of  the  shape  factor  models. 

[12]  Given  the  water-leaving  radiance  at  three  wave¬ 
lengths,  equation  (16)  still  cannot  be  solved  for  the  “big 
three”  IOPs,  ^(X,),  a^X,),  and  bbt(\ ),  because  each 
measurement  of  RSRa(\ )  yields  an  equation  with  three 
unknown  IOPs.  However,  it  is  easy  to  show  that  a  consistent 
solution  is  available  by  introducing  spectral  models  for 
aPh(K),  and  bh((\i)  [Hoge  and  Lyon ,  1996;  Hoge  et 

al .,  1999a,  1999b].  Substitution  of  such  spectral  models  for 
aPh(\)>  £*d(\),  and  bhi(\)  into  equation  (16)  yields 


+  ^(Xrf)  expf-SCX,  -  \,)1  +  MM  (Q  y 
=  -  aw(\)  -  bUK)  V-Dl-  Bf.  (17) 

Equation  (17)  now  has  only  three  unknowns,  aph(Xg)y 
and  bbt(\b)y  so  that  the  system  is  solvable  given 
measurements  of  RSRa(\)  at  three  wavelengths.  This  linear 
form  of  the  radiative  transfer  equation  remains  exact  and 
precise,  but  the  uncertainty  in  the  retrieved  IOPs  is  now 
additionally  influenced  by  the  uncertainty  in  the  IOP  models 
[Hoge  et  al.y  1999a,  1999b]  (in  addition  to  the  shape 
factors). 

[13]  At  their  respective  reference  wavelengths,  X^,  \dy  and 
\b,  the  IOPs  aph(\g),  aJ\j)y  and  bbt(\b)y  are  linearly  related 
to  the  column  matrix,  or  vector,  containing  the  hydrospheric 
constants  (sea  water  absorption  and  backscattering),  radi¬ 
ances,  and  the  shape  factors.  (It  is  easy  to  see  from 
equation  (16)  that  it  is,  in  principle,  possible  to  concurrently 
solve  for  the  radiance  derivative  term  DL(Xj)  and/or  for  the 
fractional  forward  scattering  coefficient,  Bf  in  addition  to 
the  IOPs,  given  measurements  of  RSRa(\)  at  additional 
wavelengths.  However,  the  IOP  models  then  must  be  of 
sufficient  accuracy  at  yet  a  fourth  and/or  fifth  wavelength, 


and  the  required  wavelength  dependency  of  these  models 
was  not  a  focus  of  this  study.  Therefore  such  retrievals  are 
beyond  the  scope  of  this  initial  RTE  inversion  work.) 
Equation  (17)  is  very  similar  to  the  one  used  to  analyze 
the  effect  of  radiance  errors  and  model  uncertainties  upon 
IOPs  [Hoge  and  Lyon ,  1996],  and  to  retrieve  IOPs  from 
atmospherically  corrected  airborne  and  satellite  upwelling 
radiances  [Hoge  et  al.y  1999a,  1999b,  2001]  when  retrieved 
by  semianalytic  radiance  model  inversion.  Equation  (17) 
can  therefore  be  written  in  matrix  form  as 

Dp  =  h  (18) 

Here  the  hydrospheric  vector  h  is  given  by  the  right-hand 
side  of  equation  (17).  The  IOP  state  vector  is  p  =  [aph(Kg)y 
tf^Xrf),  bbl(\b)]T,  where  T  denotes  the  transpose  and  D  is  the 
data-model  matrix  [Hoge  and  Lyony  1996;  Hoge  et  al.y 
1999a,  1999b,  2001],  which  also  contains  shape  factors. 
The  IOPs  are  immediately  determined  from  p  =  D~]h 

[14]  The  uncertainties  in  the  IOP  state  vector  p  can  be 
analyzed  in  a  manner  similar  to  other  linear  inversions 
[Hoge  and  Lyony  1996].  Since  p  =  D~]hy  both  D  and 
h  determine  p  and  the  errors  that  propagate  into  p  Because 
the  backscattering  shape  factor  fb  is  always  found  within  Dy 
fb  influences  the  propagation  of  errors  into  the  IOPs  more  so 
than  the  remaining  factors.  The  discussion  of  the  uncertain¬ 
ties  in  the  IOP  state  vector  p  caused  by  possible  singularity 
of  D~ 1  and  by  perturbations  in  D  somewhat  parallels  a 
similar  previous  discussion  [Hoge  and  Lyony  1996]  and  is 
briefly  addressed  in  a  later  section  herein. 

4.  Models  for  Shape  Factors  and 
Related  Quantities 

[15]  Inversion  of  equation  (17)  using  remotely  sensed 
ocean  color  data  requires  knowledge  of  Vy  DLy  and  Bf.  These 
quantities  in  turn  depend  on  the  shape  factors  fb  and  fLy  the 
diffuse  attenuation  coefficient  ky  and  the  mean  cosine  ratio 
R^y  as  defined  above.  For  ease  of  comparison  with  previous 
work  on  shape  factors  [Weidemann  et  al.y  1995]  and  to 
reveal  the  underlying  physics  as  much  as  possible,  explicit 
models  for  fb>fu  K  and  R ^  are  given  rather  than  DL  and  Bf. 
For  notational  convenience,  let  Xt  with  i  =  1,  2,  3,  and  4 
denote  fby  fLy  ky  and  R ^  respectively. 

[16]  Because  shape  factors,  diffuse  attenuation  functions, 
and  mean  cosines  all  depend  on  the  ambient  radiance,  they 
depend  implicitly  on  the  solar  zenith  angle  and  viewing 
direction,  as  well  as  on  the  IOPs.  The  solar  angle  and 
viewing  direction  are  known  in  any  particular  remote 
sensing  situation;  these  geometric  quantities  are  thus  avail¬ 
able  for  modeling  the  X{  in  terms  of  known  quantities. 
However,  the  IOPs  are  unknown.  Explicit  inversions  of  the 
RTE  (to  obtain  the  IOPs)  excludes  the  IOPs  from  the  models 
for  the  Xt.  However,  an  implicit,  or  iterative,  inversion  of  the 
RTE,  can  include  the  retrieved  IOPs  in  the  Xt  models,  for  the 
following  reason.  In  an  iterative  inversion,  one  starts  with 
an  initial  guess  for  the  Xiy  derived  either  from  models  that  do 
not  include  IOPs  or  from  physical  intuition.  (For  example,  a 
reasonable  initial  guess  for  fb  would  be  1,  the  value 
corresponding  to  a  constant  phase  function.  Similarly, 
fL  =  1,  k  =  0,  and  R^  =  1  would  be  acceptable  initial 
guesses.)  Using  the  initial  guesses  for  the  Xiy  the  RTE  is 
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Table  1.  Parameters  and  Their  Values  Used  to  Generate  the  Original  Database  of  184,800  Records® 


Parameter 

Values  Used  in  Hydrohght  Runs 

chlorophyll  concentration,  Chi 

0. 

1,  0.5,  1.0,  5.0,  10.0  mg  m‘3 

solar  zenith  angle,  0, 

0.1 

0,  10,  20,  30,  40,  60  degrees 

wind  speed 

0, 

10  m  s_1 

cloud  cover 

0, 

100%,  i.e.,  clear  sky  and  solid  overcast 

wavelength,  X 

412,  426,  440,  465,  490,  522.5,  555,  612,  670,  685  nm 

polar  viewing  angle  0V  (in  water,  relative  to  the  zenith) 

0.1 

D,  10,  20,  30,  40,  50,  60  degrees 

azimuthal  viewing  angle  4>v  (relative  to  Sun) 

0, 

90,  1 80  degrees 

depth,  z 

o, 

0.5,  1,  1.5,  2,  2.5,3,  3.5,4,  4.5.5  m 

"The  values  shown  in  bold  correspond  to  the  1,500  records  in  the  remote  sensing  database. 


inverted  to  obtain  initial  values  for  aph(\g),  cidO'd),  and 
bbt(\[j)y  from  which  the  IOPs  a  =  at  +  aw  =  aph  +  ad  +  aw 
and  bb  =  bbt  +  bbw  can  be  obtained  at  all  wavelengths  via 
the  IOP  models  seen  in  equation  (17).  The  Xt  models 
developed  below  are  based  on  an  assumed  phase  function 
for  particle  scattering.  Taking  the  particle  phase  function  as 
known,  the  total  constituent  (particle)  backscattering  frac¬ 
tion  Bt  =  bbt/bt  is  also  known.  Thus  the  total  constituent 
forward  scattering  coefficient  bp  can  be  obtained  from  the 
recovered  bbt  and  Bt :  bfi  =  bbt(  1  -  Bt).  The  total  beam 
attenuation  coefficient  is  then  known  from  c  =  a  +  bj)  +  bbt  + 
bw.  Therefore  models  are  developed  for  the  Xt  that  depend 
both  on  the  known  geometrical  (viewing  direction  (sub¬ 
script  v),  solar  direction  (subscript  s),  and  physical  (wind 
speed,  wavelength))  parameters,  as  well  as  on  certain  IOPs 
(namely,  a ,  c,  and  bb). 

4.1.  Database 

[n]  To  begin  the  analysis,  120  Hydrolight  [Mobley, 
2001a,  2001b]  runs  were  made  using  its  IOP  model  for 
case  1  waters  and  the  Petzold  “average  particle”  phase 
function  [Mobley  et  al. ,  1993]  for  scattering  by  the  particles. 
This  case  1  IOP  model  is  a  two-component  model,  pure 
water  plus  “everything  else.”  The  non-water  absorption  and 
scattering  coefficients  are  parameterized  in  terms  of  the 
chlorophyll  concentration  according  to  commonly  used 
models  by  Mobley  [1994,  equations  (3.27)  and  (3.40)]. 
The  input  for  these  runs  covered  a  wide  range  of  chloro¬ 
phyll  concentrations,  solar  zenith  angles,  cloud  covers,  and 
wind  speeds.  Each  Hydrolight  run  generated  output  at 
various  wavelengths,  depths,  and  viewing  directions.  The 
resulting  database  potentially  contains  millions  of  records, 
where  one  record  corresponds  to  a  particular  set  of  input 
values,  output  values  for  a  particular  viewing  geometry, 
wavelength,  depth,  etc.,  and  the  values  of  the  four  Xt.  Some 
of  these  records  are  not  of  great  interest,  for  example, 
records  whose  azimuthal  viewing  directions  <j>v  differ  by 
only  1 5  degrees  (the  resolution  of  <J>V  in  the  standard  version 
of  Hydrolight).  Therefore  selected  records  were  used  to 
generate  a  database  of  more  manageable  size  but  one  that 
still  covers  the  range  of  parameter  values  relevant  to  most 
remote  sensing.  Table  1  shows  the  input  and  output  values 
in  this  database,  which  was  used  in  the  initial  investigation 
of  the  functional  forms  of  the  Xt.  Each  of  the  four  Xt  was 
computed  for  each  parameter  combination  represented  in 
Table  1. 

[is]  First,  the  sensitivity  of  the  Xt  to  the  various  param¬ 
eters  (wind  speed,  viewing  direction,  IOPs,  etc.)  available 
for  construction  of  models  for  the  X{  was  examined.  It  was 


found  that  the  surface  wind  speed  has  a  negligible  effect  on 
each  of  the  Xt  (less  than  1  %  difference  in  Xt  for  the  0  and 
10  m  s"1  wind  speeds,  with  all  else  being  equal).  Thus  the 
wind  speed  was  not  considered  in  subsequent  modeling. 
Likewise,  for  similar  reasons,  only  the  clear-sky  data,  which 
are  of  greatest  interest  for  remote  sensing  applications  were 
included.  As  noted  above,  for  remote  sensing  applications 
the  shape  factors  Xt  need  evaluation  only  at  depth  z  =  0. 
(Note,  however,  that  the  values  of  the  Xt  at  z  =  0  incorporate 
the  effects  of  all  the  absorption  and  multiple  scattering 
occurring  throughout  the  entire  water  column.)  Thus  only 
output  from  the  Hydrohght  runs  at  z  =  0  was  retained.  The 
original  database  included  records  generated  for  azimuthal 
viewing  angles  of  4>v  =  0  (looking  toward  the  Sun)  and  <J>V  = 
180  degrees  (looking  away  from  the  Sun).  Most  remote 
sensing  is,  or  can  be,  accommodated  at  azimuthal  angles  of 
4>v  «  90  degrees,  which  minimizes  Sun  glint  and  instrument 
self  shading.  Thus  only  the  records  corresponding  to  <j>v  = 
90  degrees  were  included.  Likewise,  remote  sensing  gener¬ 
ally  uses  in-air  nadir  viewing  directions  0Vfl  of  less  than 
60  degrees,  which  correspond  to  in-water  angles  0V  < 
40  degrees.  Eliminating  the  larger  in-water  nadir  viewing 
angles  (0V  =  50  and  60  degrees  in  Table  1)  gives  a  final 
“remote  sensing”  data  set  of  1500  records,  which  was  used 
to  determine  models  for  the  Xt.  The  parameter  values 
corresponding  to  this  remote  sensing  database  are  shown 
in  bold  in  Table  1 . 

4.2.  Determining  Functional  Forms 

[19]  Let  Pk  with  k  =  1,...,A^,  denote  the  parameters 
(wavelength,  viewing  direction,  IOPs,  etc.)  to  be  used  in 
modeling  the  Xt.  These  parameters  include  those  seen  in 
Table  1 ,  as  well  as  the  absorption,  scattering,  and  backscat¬ 
tering  coefficients  (which  are  functions  of  the  chlorophyll 
concentration  if  case  1  water  is  assumed). 

[20]  The  simplest  possible  model  for  the  Xt  is  a  linear 
function  of  the  Pk: 

Nk 

X  =  Y,  <*kPk,  «=  1.....4,  (19) 

*=I 

The  <y.ik  are  fitting  coefficients  whose  values  are  to  be 
determined;  a  different  set  of  coefficients  is  needed  for  each 
factor  Xt.  That  is,  a  large  linear  least  squares  problem  was 
initiated  to  determine  if  this  model  was  adequate  to  fit  the 
various  factors.  Not  surprisingly,  the  fits  were  unsatisfac¬ 
tory.  In  other  words,  the  ocean  is  more  complicated  than 
equation  (19). 
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Figure  1.  Sun,  sensor,  and  single-scattering  geometry.  Equation  (1)  measures  the  polar  angle  0  from  the 
nadir  with  +z  downward,  and  Hydrolight  computes  in  this  coordinate  system  as  well,  but  the  above¬ 
surface  solar  and  viewing  directions  are  shown  for  convenience  as  being  measured  from  the  zenith. 


[21]  The  goal  then  became  to  develop  models  that  are 
nonlinear  in  the  parameters  but  still  reflect  the  underlying 
physics  of  the  Xt.  Thus  replace  equation  (19)  by 

N(.) 

=  W»-  (20) 

This  notation  means  that  each  function  Fik  will  contain 
some  subset  of  the  parameters  Pk  and  will  have  its  own  set 
of  fitting  coefficients.  For  example,  as  shall  be  seen  in  the 
next  section,  the  model  for  fb  will  contain  terms  involving 
the  scattering-to-backscattering  ratio  and  a  nonlinear 
function  of  the  three  geometric  parameters  (05,  0V,  and  4>v), 
with  four  fitting  coefficients  in  all. 

4.3.  Model  for  fb 

[22]  To  obtain  initial  guidance  about  the  possible  func¬ 
tional  form  for  an  fh  model,  single-scattering  theory  was 
first  used  to  evaluate  the  definition  of  fb  seen  in 
equation  (3).  According  to  the  single-scattering  approxi¬ 
mation  (SSA),  the  downwelling  radiance  is  [ Gordon , 
1994,  equation  (1.30)] 

-ct/m, 

> 


where  p,  =  cos  0  and  (p.v,  4^)  denote  the  in-water  direction  of 
the  Sun’s  direct  beam  as  refracted  through  a  level  sea 
surface;  6  is  the  Dirac  delta  function.  Inserting  equation  (21) 
into  equation  (3),  integrating,  and  evaluating  the  result  at 
z  =  0  gives 

fb  «  (22) 

Note  that  for  isotropic  scattering,  (3  =  1/4tt,  b  =  2 bbi  giving 
fb  =  1,  as  expected.  As  the  phase  function  becomes  more 
anisotropic,  b/bb  increases,  and  fb  increases.  In  equation  (22) 
the  phase  function  indicates  that  to  first  order,  fb  involves 
downwelling  radiance  that  is  singly  scattered  from  the  Sun’s 
direct  beam  into  the  viewing  direction.  In  remote  sensing 
applications  the  total  phase  function  (water  plus  particles)  will 
be  unknown,  but  for  any  given  phase  function  the  contribu¬ 
tion  to fb  will  depend  on  the  scattering  angle  vj/  corresponding 
to  the  Sun’s  downward  beam  being  scattered  into  the 
upward  viewing  direction.  This  scattering  angle  is 
equivalent  to  the  easily  computed  Sun  sensor-included 
angle  6,,  as  shown  in  Figure  1 .  Given  the  Sun’s  location  (05; 
<$>s  =  0),  the  viewing  direction  (0V,  4>v),  £  is  given  by 

cos  f;  =  cos  0S  cos  0V  +  sin  0S  sin  Gy  cos  <|>v . 


4 >\z)  =  Eo, 6(p'  -  m.,)6(4>'  -  4>> 


(21) 
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1.30 


□  b/bb  £  30 
A  30  <  b/bb  ^  45  J 
X  45  <  b/bb 


A  1  -*  4.«.  111.1  11  J  i  I  |  , 


0  10  20  30  40  50  60  70  80 

f  [deg] 

Figure  2.  Values  of  fb  plotted  as  a  function  of  Sun  sensor  included  angle  6,  and  coded  to  show  the 
dependence  on  blbh. 


Thus,  as  a  preliminary  functional  form  for  the  fb  model, 
consider  the  estimated  backscattering  shape  factor 

fb  «  ^=(0- 


where  S(£)  is  a  function  of  angle  £  whose  form  is  to  be 
determined. 

[23]  Figure  2  shows  the  1,500  values  of fb  in  the  remote 
sensing  database  plotted  as  a  function  of  £  and  b/bb.  This 
figure  suggests  that  a  cosine  function  may  capture  the  £ 
dependence.  Thus  construct  a  model  of  the  form 

fb  =  a*  +  ot2  ^  cos(a3  0,  (23) 

where  oq,  a2,  and  a3  are  fitting  coefficients  (the  a*  of 
equation  (20)  for  model  i  =  1)  whose  values  are  to  be 
determined  by  minimizing  the  squared  difference  between 
fb  and  fh  for  the  1,500  values  in  the  remote  sensing 
database.  Note  from  the  points  in  Figure  2  that  fb  is  largest 
for  small  b/bbi  and  vice  versa,  which  is  contrary  to  the 
behavior  predicted  by  equation  (22).  This  reversal  may  be 
due  to  the  dominance  of  multiple  scattering  in  ocean  waters, 
but  further  investigations  would  be  necessary  to  understand 
this  discrepancy  between  the  SSA  predictions  and  the 
Hydrolight  predictions,  which  include  all  orders  of  multiple 
scattering  and  other  effects  not  included  in  the  SSA.  In  any 
case,  there  is  a  clear  dependence  on  b/bhi  which  can  be 
modeled. 

[24]  The  best  fit  coefficients  otj  in  equation  (23)  were 
determined  by  least-squares  minimization  using  a  variety  of 
numerical  techniques  appropriate  for  nonlinear  functions. 
After  comparing  the  model  predictions  fb  with  the  actual 
fb  values,  it  was  seen  that  not  all  of  the  blbb  dependence  was 


captured  by  the  model  of  equation  (23).  Some  experimen¬ 
tation  showed  that  the  remaining  blbh  dependence  could  be 
accounted  for  by  adding  another  term  proportional  to  b/bb. 
Thus  the  final  fb  took  the  form 


fb  =  on  +  012^-  cos(ot3  O+cu  £ 


=  01+04  — 

Ob 


1  +  —  cos(o3  0 
04 


(24) 


Equation  (24)  shows  that  the  additive  term  is  equivalent  to 
keeping  the  general  form  of  the  model  suggested  by  the 
SSA,  but  picking  a  different  angular  function  Eft).  The 
complicated  dependence  of fb  on  b/bb  and  £  is  not  surprising 
if  one  remembers  that  the  S(0  function  is  fundamentally  an 
attempt  to  parameterize  the  unknown  phase  function  effects 
for  a  given  Sun  and  viewing  geometry;  thus  the  scattering 
and  geometric  effects  are  not  independent.  The  final  set  of 
fitting  coefficients  for  the  model  in  equation  (24)  is  shown 
in  Table  2. 

[25]  Figure  3  shows  the  model  and  actual  fh  values.  The 
dashed  lines  are  the  5%  error  bounds.  Using  the  model  of 
equation  (24),  96.3%  of  the  predicted  values  are  within  5% 
of  the  correct  value;  the  linear  correlation  coefficient 
between  the  model  and  actual  points  is  r  =  0.955.  There 
is  no  systematic  dependence  on  b/bb  or  6,  of  the  model 


Table  2.  Best  Fit  Coefficients  for  the  fh  Model  of  Equation  (24) 


Coefficient 

Value 

a  1 

1.2077 

012 

0.001977 

015 

3.3790 

<*4 

-0.004863 
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Figure  3.  Comparison  of  modeled  and  actual  fb  values,  using  the  model  of  equation  (24).  Here  96.3%  of 
the  modeled  values  lie  within  the  dashed  lines,  which  represent  values  with  5%  of  the  correct  value;  the 
model-actual  correlation  coefficient  is  r  -  0.955. 


discrepancies  seen  in  the  individual  points  of  Figure  3.  It 
would  be  possible  to  continue  adding  ad  hoc  terms  in  other 
variables  to  equation  (24)  and  perhaps  reduce  the  model- 
data  discrepancy  even  more.  However,  such  a  process  is 
likely  to  deviate  from  physical  foundations,  with  the  end 
result  that  the  final  model  would  not  be  applicable  beyond 
the  exact  conditions  used  to  generate  the  present  remote 
sensing  data  set.  For  this  initial  study  it  is  best  to  be  content 
with  the  simple  model  of  equation  (24). 

4.4.  Model  for  bf[  1  —  fL\ 

[26]  Since  bf  is  an  IOP  the  model  required  for  Bf\sfL.  Just 
as  with  fbi  use  the  SSA  for  guidance  as  to  the  general  form 
of  the  model  for  fa.  In  the  SSA  the  upwelling  radiance  is 
[ Gordon ,  1994,  equation  (1.32)] 

=  -£rf(0)P(M,1+J,|i'1  ♦')— —e-a^.  (25) 

c  M-j  M* 


function  was  parameterized  in  terms  of  the  Sun  sensor 
included  angle  £  via  a  model  of  the  form 

Unfortunately,  plots  offi  similar  to  Figure  2  did  not  suggest 
a  clear  functional  form  for  3(0  or  show  any  significant 
dependence  on  b/c.  This  failure  of  the  SSA  to  provide  a 
functional  form  ioifL  is  not  suiprising  because  inherently 
involves  at  least  two  scatterings,  and  multiple  scattering  can 
be  expected  to  make  an  important  contribution  to  the 
upwelling  radiance. 

[27]  Linear  correlations  between^  and  the  various  avail¬ 
able  fitting  parameters  were  then  examined.  The  results  are 
seen  in  Table  3. 

[28]  The  only  potential  model  parameter  that  correlates 
well  with  fL  is  the  solar  zenith  angle  Qs.  A  plot  0  ffL  versus  0y 


Note  that  [is  >  0  and  \i  <  0.  Inserting  this  SSA  radiance 
into  equation  (4),  the  definition  of fL)  integrating,  and  setting 
z  =  0  gives 


,  &2£/(0) 
JL  ~  bfcVjo 


<t>'»  M,  Mm  4>;) ' — — ~,d\>!  dtf. 


M-j-M- 


In  most  ocean  waters,  b  «  bf.  However,  further  simplifica¬ 
tion  is  difficult.  The  remaining  integrals  describe  how  the 
Sun’s  down  welling  direct  beam  is  first  scattered  upward  and 
then  scattered  again  into  the  viewing  direction.  The  most 
that  can  be  said  is  that  this  is  some  function  of  the  scattering 
phase  function  and  the  viewing  geometry.  As  with  fb ,  this 


Table  3.  Correlation  Between  the  Radiance  Shape  Factor  fL  and 
Various  Parameters 


Parameter 

Correlation  Coef.  r  With  fL 

Absorption  coefficient  a 

0.291 

Scattering  coefficient  b 

0.205 

Albedo  of  single  scattering  b/c 

-0.076 

Backscattering  coef  bb 

0.180 

Forward  scattering  coef  bf 

0.206 

Backscattering  fraction  bjb 

-0.343 

Forward  scattering  fraction  bf/b 

0.343 

Backscattering  to  absorption  bja 

-0.231 

Wavelength  X 

0.210 

Solar  zenith  angle  0, 

0.778 

Polar  viewing  angle  0V 

-0.006 

Sun  sensor  included  angle  £ 

-0.562 
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actual  /l 

Figure  4.  Comparison  of  modeled  and  actual^  values,  using  the  model  of  equation  (24).  Here  93.0% 
of  the  modeled  values  lie  within  the  dashed  lines,  which  represent  values  with  2%  of  die  correct  value; 
the  model-actual  correlation  coefficient  is  r  =  0.829. 


suggested  a  sine  function  for  0*  (although  a  simple  linear 
function  is  almost  as  good).  Thus  a  model  of  the  form  fL  = 
a5  +  ot$  sin(a70>r)  was  tried.  The  residuals  of  this  model 
showed  a  weak  wavelength  dependence.  After  considerable 
experimentation  the  model  chosen  was 

/i  =  o5  +  sin(o7  6,) .  (26) 


errors.  The  dotted  lines  in  Figure  5  thus  show  absolute 
errors  of  ±0.02  m-1;  94.5%  of  the  Bf  have  errors  smaller 
than  this.  The  model-actual  correlation  coefficient  is  r  = 
0.966. 

4.5.  Model  for  k 

[30]  The  SSA  again  suggests  a  functional  form  for  the  £(0, 
<t>,  z  -  0)  model.  Differentiating  the  SSA  up  welling  radiance 
of  equation  (25)  gives 


The  best  fit  values  of  the  coefficients  are  a5  =  1.0247,  = 

0.4584  (for  X  in  nanometers),  and  a7  =  0.1221  (for  0* 
measured  in  radians).  Figure  4  shows  the  scatterplot  for  this 
model.  Here  93.0%  of  the  model  predictions  are  within  2% 
of  correct  (points  lying  between  the  dashed  lines);  the 
correlation  coefficient  is  r  -  0.829.  Although  it  is  possible 
to  obtain  slightly  better  fits  by  including  IOPs  in  an  ad  hoc 
fashion,  the  model  of  equation  (26)  was  selected  because  of 
its  simplicity  and  because  of  the  lack  of  physical  guidance 
for  the  IOP  dependence. 

[29]  Although  one  is  unable  to  model  the  remaining 
variability  of  fL  in  terms  of  the  IOPs  or  other  parameters, 
this  may  be  of  little  importance  in  predicting  fL  itself 
because  fL  is  always  near  1 .  Perhaps  more  important  is  the 
fact  that  the  variability  in  fL  determines  the  variability  in  the 
fractional  forward  scattering  coefficient  Bf  -  bf  (1  —  fL). 
Small  fractional  errors  in  fL  can  cause  large  fractional  errors 
in  Bf.  Figure  5  shows  the  resulting  scatterplot  for  Bf, 
computed  using  the  exact  values  of  bf  as  found  in  the 
database.  Although  93.0%  of  the^  values  are  within  2% 
of  their  correct  value,  only  5.4%  of  the  Rvalues  are  within 
2%  of  the  correct  value;  58.1%  of  the  Bf  are  within  20%  of 
the  correct  value.  However,  a  percentage  error  criterion  may 
be  misleading  for  Bf  because  of  the  cluster  of  points  near 
zero,  where  small  absolute  errors  can  be  large  fractional 


*(«,♦,  0) 


1  dLu(z)  1  c 

Lu(z)  dz  JI=0  cos(0,) 


Although  k  is  strongly  correlated  with  a  (r  =  0.99),  it  is  not 
strongly  correlated  with  c  (r  =  0.56).  At  the  level  of  the 
quasi-single-scattering  approximation  (QSSA),  which  often 
works  well  for  upwellmg  radiances,  c  «  a  +  bbi  in  which 
case,  k  «  (a  +  bb)/ cos(0v).  This  suggested  a  model  k  with  the 
functional  form 


-  _  as  a  +  (*9  bb 
cos(0j) 


(27) 


When  equation  (27)  was  used  to  fit  the  points  in  the  remote 
sensing  data  set,  the  points  separated  into  distinct  groups  for 
0,  >  60°  and  for  0*  <  60°.  Plots  of  k  and  the  residual  error  in 
k  as  functions  of  05  suggested  that  an  additive  sin(0.y)  term 
would  represent  the  0V  dependence  better  than  the  l/cos(05) 
factor  in  equation  (27).  This  then  gave  the  final  k  model: 


k  =  cx%  a  +  ol$  bi,  +  aio  sin^) .  (28) 

The  best  fit  parameters  are  a8  =  1 .0896,  a9  =  —0.593 1,  and 
aio  =  0.0492  (for  a ,  bbi  and  k  in  inverse  meters).  (Note  that 
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Figure  S.  Comparison  of  modeled  and  actual  Rvalues,  computed  using  th  zfL  model  of  equation  (26). 
The  dashed  lines  are  20%  error  bounds;  the  dotted  lines  are  ±0.02  m-1  error  bounds. 


although  0V  in  equation  (28)  refers  to  the  polar  angle  of  the 
solar  beam  in  water,  the  in-air  solar  zenith  angle  for  0,,  can  be 
used  because  the  index-of-refraction  factor  that  converts 
sin^  in  air)  to  sin(05  in  water)  is  incorporated  into  otio)  The 
points  still  separate  somewhat  by  0*  >  60°  and  for  0j  <  60° 
but  not  as  much  as  for  equation  (27).  Although  equation  (28) 
has  lost  some  of  its  intuitive,  first-order  physics,  namely  the 
l/cos(0j)  factor  in  equation  (27),  the  final  model  does  a  better 
job  of  predicting  k ,  which,  of  course,  is  influenced  by 
multiple  scattering  and  other  effects  not  included  in  the  SSA. 


Figure  6  shows  the  scatterplot  for  the  k  model.  Because  of 
small  differences  near  k  =  0,  only  76.8%  of  the  points  are 
within  20%  of  the  correct  value.  However,  93.6%  of  the 
points  lie  within  an  absolute  error  of  ±0.05  m-1.  The  model- 
actual  correlation  coefficient  is  r  -  0.993. 

4.6.  Model  for  R ^ 

[31]  A  model  for  =  jj^in  ai^/fl^in  water  at  z  =  0)  can 
be  constructed  simply  by  using  Snell's  law  to  refract  the 
direct  solar  beam  through  a  level  water  surface.  The  result  is 


0.0  0.2  0.4  0.6  0.8 

actual  k  [m-1] 


Figure  6.  The  k  model  of  equation  (28).  The  dashed  lines  are  the  20%  error  bounds,  and  the  dotted  lines 
are  ±0.05  m-1  error  bounds. 
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Figure  7.  The  model  of  equation  (29)  applied  to  the  RS  data  set.  The  dashed  lines  are  5%  error  bounds. 


K  =  an- 


COS  0y 


/  sin  0,\ 1 

cos 

sin  1 

{—)\ 

=  an0(0J), 


(29) 


where  n  =  1.34  is  the  index  of  refraction  of  water.  The 
ratio  R ^  is  independent  of  the  viewing  direction  because 
the  mean  cosines  are  computed  from  integrals  over  the 
direction  of  the  radiance  distribution;  there  are  conse¬ 
quently,  many  fewer  distinct  points  in  the  data  set.  When 
equation  (29)  is  applied  to  the  RS  data  set,  the  best  fit 
value  is  an  =  0.869.  Figure  7  shows  the  results  of  this 
model  applied  to  all  points  in  the  RS  data  set.  The  six 
groups  of  points  correspond  to  the  six  solar  angles  in  the 
data  set:  0V  =  0,  10,  20,  30,  40,  and  60  degrees.  Here 
85.7%  of  the  points  lie  within  5%  of  the  correct  value,  as 
shown  by  the  dashed  lines.  The  correlation  coefficient  is 
r  =  0.964. 


5.  Remotely  Sensed  Reflectance  RSRa  Estimation 

[32]  Models  that  estimate  the  shape  factors  and  related 
quantities  with  varying  degrees  of  certainty  are  now  com¬ 
plete.  The  question  next  arises  as  to  how  well  one  can 
predict  the  remotely  sensed  reflectance  RSRa  using  these 
models  within  equation  (9).  Figure  8  gives  the  answer: 
68  . 1%  of  the  RSRa  predictions  fall  within  10%  of  the  correct 
values,  and  95.7%  fall  within  20%  of  correct  (shown  by  the 
dashed  lines  in  Figure  8);  85.6%  of  the  predictions  are 
within  ±0.0005  sr_1  of  the  correct  value.  The  model-actual 
correlation  coefficient  is  r-  0.983. 


6.  Discussion 

[33]  To  facilitate  a  brief  comparison  of  the  shape  factor 
models,  propagation  of  errors  into  the  retrieved  IOPs,  and 


discussion  of  future  inversion  research,  all  the  models  are 
reassembled  below. 


A 


=  ai  +  04 


b_ 

h 


1  -I-  —  cos(a3  £) 
04 


(30) 


Bf  =  bf[\  - }l ) 


where 


h 


=  a5  +  a« 


sm(ot70j).  (31) 


k 


=  ag  a  +  cu)  bh  4-  otio  sin(0j). 


(32) 


K  =  an 


cos  Oj 

cos 

sin-1  | 

(*?) 

4 

0116(65). 


(33) 


Although  they  were  derived  from  a  physical  basis,  it  was 
seen  that  the  models  could  take  various  forms.  At  this  early 
stage  of  development  the  above  models  probably  represent 
the  starting  point  of  their  eventual  evolution. 

[34]  The  highly  important  fb  model  contains  (1)  two 
IOPs:  bb  and  bf  (but  in  a  ratio  combination  b/bb  =  [(&*,  + 
b/)lbb]  =  [1  +  bf/bb]),  (2)  the  most  model  coefficients  (four), 
and  (3)  the  Sun  sensor  included  angle  6,  (but  not  the  solar 
zenith  angle  0*  as  do  all  the  other  models).  In  contrast,  the  Bf 
model  contains  (1)  the  solar  zenith  angle  and  one  IOP  (bj) 
and  (2)  the  sole  wavelength  dependence  found  within  the 
models.  The  k  model  contains  only  one  IOP,  bb,  and  the 
solar  zenith  angle,  05.  The  model  contains  no  IOPs;  only 
the  solar  zenith  angle  0V.  (Inversion  of  the  shape  factor  RTE 
also  requires  models  for  those  IOPs  that  are  to  be  retrieved. 
For  example,  the  phytoplankton  absorption  coefficient  aph , 
the  CDOM/detritus  absorption  coefficient  ad  and  total 
constituent  backscattering  bbt  as  given  in  equation  (17). 
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0.000  0.002  0.004  0.006  0.008  0.010 
actual  RSRa  [sr-1] 

Figure  8.  RSRa  as  modeled  by  equation  (9)  using  the  four  shape  factor  models.  The  dashed  lines  are 
20%  error  bounds;  the  dotted  lines  are  ±0.0005  sr-1. 


These  IOP  models  [Hoge  and  Lyon ,  1996]  are  considered 
more  mature  than  the  shape  factor  models.  Uncertainty 
propagated  into  retrieved  IOPs  by  the  IOP  models  used 
within  a  semianalytic  radiance  model  inversion  has  been 
studied  [Hoge  and  Lyon ,  1996].) 

[35]  Thus,  to  initiate  an  iterative  inversion,  starting  values 
are  required  for  both  b b  and  bf.  Physics  demands  that  bb  > 
bhwy  where  bbw  is  the  backs cattering  coefficient  for  water. 
One  possible  method  for  selecting  the  starting  value  for  bb  is 
to  retrieve  it  by  first  executing  a  semianalytic  model 
inversion  [Hoge  and  Lyon ,  1996;  Hoge  et  al. ,  1999a, 
1999b,  2001].  Then  it  can  continually  be  updated  after  each 
shape  factor  RTE  inversion  in  equation  (17)  since  bb  -  bbw  + 
bbt.  Similarly,  physics  dictates  and  limits  the  range  of  bf  fox 
the  first  iteration  of  the  shape  factor  RTE  inversion:  b/>  b ^ 
where  bf w  is  the  forward  scattering  coefficient  for  water. 
Although  bf=  bfr  can  perhaps  be  used  as  the  starting  value 
for  the  first  iteration,  future  research  efforts  must  develop 
methods  for  better  (1)  selection  of  starting  values  and  (2) 
updating  of  the  value  during  subsequent  iterations.  Like  bby 
the  bf  can,  in  principle,  be  retrieved  using  equation  (17). 
This  too,  however,  presents  some  concerns:  (1)  few  if  any 
models  exist  for  bf  to  allow  its  retrieval  by  equations  (17) 
and  (2)  a  concurrent  retrieval  of  bf  potentially  weakens  the 
retrieval  of  the  desired  aphi  ad ,  and  bbt.  Detailed  error 
propagation  analyses  of  the  shape  factor  RTE  inversion 
are  outside  the  scope  of  this  present  paper,  but  a  brief 
discussion  of  the  relative  influence  of  the  shape  factor 
models  on  the  desired  IOP  state  vector,  p  =  [aphQ^g),  a<*(\/) 
bhtQ'b)]7*  is  provided  in  the  following  section. 

6.1.  Uncertainties  in  the  IOP  State  Vector  p 

6.1.1.  Sensitivity  of  p  to  Perturbations  in  the 
Data-Model  Matrix  D 

[36]  As  already  noted,  the  inversion  of  the  shape  factor 
form  of  the  RTE  is  exact  from  the  standpoint  of  radiative 


transfer  theory,  and  uncertainties  in  the  retrieved  IOPs 
within  the  IOP  state  vector  p  are  due  only  to  the  accuracy 
of  the  (1)  shape  factor  models  and  their  related  quantities 
and  (2)  IOP  models.  Perturbations  within  D  arise,  for 
example,  from  the  water-leaving  radiances,  scalar  irradian- 
ces,  IOP  models,  and  backscattering  shape  factor  contained 
within  it.  Similarly,  uncertainties  in  h  arise  from  the 
radiances,  irradiances,  hydrospheric  constants  (or  IOP  con¬ 
stants  aw  and  bbw )  for  sea  water,  as  well  as  fb ,  dLu(\t)ldzy 
^/(\)>  /L(\)>  and  cos  0.  Relative  to  h,  the  data-model 
matrix,  Z>,  plays  the  major  role  in  the  propagation  of  errors 
»nto  p  since  \\p  - p'\\l\\p\\  <  k .(D)  (||Ad||/||D||  +  ||6*||/||6|h, 
where  ||Z>||  is  the  determinate  of  D  and  k,(Z>)  =  ||Z)||  ||Z)  7|| 
[< Ortega ,  1990;  Hoge  and  Lyon ,  1996].  The  latter  expression 
is  the  condition  number  of  />,  and  A/>  and  hh  represent 
uncertainty  or  perturbation  of  D  and  h,  respectively.  Here  p 
is  the  perturbed  solution  of  p.  The  first  expression  simply 
states  that  to  first  order  the  relative  error  in  p  can  be  k ,(/)) 
times  the  relative  error  in  D  and  h  Thus  the  propagation 
into  p  of  the  relative  errors  of  both  D  and  h  is  governed  by 
the  condition  number  of  D  For  any  norm,  1  <  *,(£>)  <  00. 
For  the  limiting  cases:  k(D)  =  1,  D  is  said  to  be  perfectly 
conditioned,  while  for  k .(D)  =  00,  D  is  singular.  For 
intermediate  values  of  k(D)  the  interpretation  of  the  condi¬ 
tion  number  is  very  subjective  and  must  be  evaluated 
separately.  For  large  k ,(D)  the  D  matrix  is  said  to  be  ill 
conditioned  and  large  errors  may  be  found  in  p  For  small 
k .(D)  the  D  matrix  is  said  to  be  well-conditioned  and  smaller 
errors  may  be  found  in  p.  Of  the  shape  factor  components 
only  fb  occurs  in  D  (via  V)  and  therefore  provides  the 
strongest  influence  on  the  IOP  retrieval  errors.  This  is  in 
agreement  with  Zaneveld  [1995],  who  concluded  that^  is 
most  critical  since  the  in-water  remotely  sensed  reflectance 
(see  equation  (2))  is  directly  proportional  to  it.  It  is  for  this 
reason  that  shape  factor  RTE  component  model  develop¬ 
ments  should  probably  focus  on  fb. 
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[37]  However,  if  other  RTE  components  such  as  DL  (or 

dLJdz)  and/or  Bf  are  solved  for,  then  they  too  will  appear  in 
the  D  matrix  and  thereby  further  increase  the  errors  in  p 
Also,  in  general,  the  condition  number  increases  as  the 
number  of  unknowns  increases  [ McCormick ,  1992],  con¬ 
tributing  still  more  uncertainty  in  p  In  part,  the  additional 
uncertainty  in  p  will  then  be  due  to  DL  and/or  Bf  model 
errors.  If  DL  and  Z^are  both  zero,  then  equation  (13)  shows 
that  both  Rrs  and  RSRa  are  linearly  proportional  to  b\J 
(< a  +  bb),  which  is  a  well  known  approximate  functional 
form  for  the  dependence  of  water-leaving  reflectances  on 
the  absorption  and  backscattering  coefficients.  However, 
other  work  [Gordon  et  al. ,  1988]  suggests,  but  does  not 
prove,  that  since  the  shape  factor  RTE  is  exact,  both  DL  and 
Bf  cannot  concurrently  be  zero.  For  example,  the  semi- 
analytic  model  contains  a  quadratic  term  [bj(a  +  bf,)]2,  and 
this  suggests  that  the  multiplier  R^  in  the  numer¬ 

ator  of  equation  (13)  above  must  jointly  account  for  both 
linear  and  quadratic  variability  in  b\J{a  +  bb)  if  both  Di  and 
Bj  are  null.  This  comparison  to  the  semianalytic  model 
[Gordon  et  al.,  1988]  also  suggests,  but  does  not  prove,  that 
(1)  Dl  and  Bf  (when  not  being  solved  for)  jointly  contribute 
only  a  small  amount  to  the  reflectances  and  thus  to  the  IOP 
retrievals  and,  (2)  accordingly,  their  contribution  to  re¬ 
trieved  IOP  uncertainty  may  not  be  strong. 

6.1.2.  Sensitivity  of  p  to  Perturbations  in  h 

[38]  While  the  condition  of  D  is  most  important  in  deter¬ 
mining  the  errors  in  the  IOP  state  vector  p ,  the  error 
propagation  equation,  || p  -  p'\\/\\p\\  <  k(D)  (||  Az>||/||/)||  + 
||6*||/||A||),  shows  that  uncertainties  in  A  also  propagate  into/?. 

6.2.  Future  Studies 

[39]  To  fully  understand  how  shape  factor  model  errors 
affect  the  accuracy  of  retrieved  IOPs,  it  is  necessary  to 
perform  in-depth  studies  of  the  iterative  shape  factor  inver¬ 
sion  algorithm  outlined  above  by  (1)  its  application  to 
synthetic  and  (2)  actual  data  sets.  Thus  future  research  in 
our  laboratories  will  study  the  details  of  the  shape  factor 
inversion  in  a  controlled  environment  such  as  Hydrolight¬ 
generated  synthetic  RSRa  data,  for  which  the  correct  IOP 
values  are  known  from  the  input  to  the  Hydrolight  computer 
program.  The  inversion  will  then  be  applied  to  actual  RSRa 
or  Rrs  field  data  that  further  contain  experimental  errors. 
Too,  the  convergence  and  condition  (i.e.,  well  conditioned 
or  ill  conditioned)  of  successive  iterations  using  the  re¬ 
trieved/updated  bf  and  bb  must  also  be  assessed.  Results  of 
these  anticipated  studies  of  the  inversion  of  the  shape  factor 
RTE  are  the  subject  of  future  publications. 

Notation 

a  total  absorption  coefficient,  at  +  aw,  m-1;  denotes 
“in  air”  when  used  as  a  subscript. 
ad  absorption  coefficient  of  CDOM  and  detritus,  m“ 1 . 
aph  absorption  coefficient  of  phytoplankton,  m”1. 
at  total  constituent  absorption  coefficient,  at  =  aph  + 
adi  m-1. 

aw  absorption  coefficient  of  water,  m"1. 
b  total  scattering  coefficient,  m-1. 
bb  total  backscattering  coefficient,  bb  =  bbw  +  bbt>  m“ 1 . 
bf  total  forward  scattering  coefficient,  m-1. 
bf)  total  constituent  forward  scattering  coefficient,  m~ 1 . 


bbt  total  constituent  backscattering  (TCB)  coefficient, 
m”1. 

bb„  backscattering  coefficient  of  seawater,  m~ 1 . 

Bf  fractional  forward  scattering  coefficient,  defined 
by  equation  (11),  m_I. 

Bt  bbt/bt,  the  total  constituent  (particle)  backscattering 
fraction. 

c  beam  attenuation  coefficient,  c  =  a  +  A,  m“!. 

CDOM  chromophoric  dissolved  organic  matter. 

D  data  and  model  matrix. 

||Z>||  determinate  of  D 

Dl  radiance  derivative  term,  defined  by  equation  (10), 
m-1. 

Eod(z)  downwelling  scalar  irradiance,  in  water,  W  m-2 
nm"1. 

Eoda  downwelling  scalar  irradiance,  in  air,  just  above 
sea  surface,  W  m-2  nm-1. 

E^z)  downwelling  plane  irradiance,  in  water,  W  m-2 
nm"1. 

Eda  downwelling  plane  irradiance,  in  air,  W  m-2 
nm-1. 

Eu(z)  upwelling  plane  irradiance,  in  water,  W  m“2  nm  - 1 . 

Eua  upwelling  plane  irradiance,  in  air,  W  m-2  nm-1. 

Fik  modeling  function;  contains  some  subset  of  the 
parameters  Pk. 

fb  backscattering  shape  factor,  dimensionless. 

fb  estimated  backscattering  shape  factor,  dimension¬ 
less. 

fL  radiance  shape  factor,  dimensionless. 

estimated  radiance  shape  factor,  dimensionless. 

fi  average  of  fL  values  having  0*  in  remote  sensing 
data  set. 

g  phytoplankton  Gaussian  model  spectral  width 
parameter,  nm. 

h  vector  of  hydrosphenc  constants,  shape  factors, 
radiance  attenuation  coefficient,  m-1. 

IOP  inherent  optical  property. 

KLu  diffuse  attenuation  coefficient  for  upwelling 
radiance,  KLv  =  -d[logZ,M(z,  0,  <j>)]/dz. 

Lu  upwelling  radiance,  below  sea  surface,  W  m-2  sr-1 
nm-1. 

Eua  upwelling  radiance,  in  air,  just  above  sea  surface, 
W  m-2  sr-1  nm-1. 

M  [(fr)/[(l  —  rR)ni]  «  0.54  for  nominal  sea 
conditions,  dimensionless. 

n  total  constituent  backscattering  coefficient  spectral 
model  exponent,  as  used  in  equation  (17), 
dimensionless. 

n„  index  of  refraction  of  sea  water,  dimensionless. 

Pk  modeling  parameter;  k  =  1,...,  Nk ,  denote  the 
parameters  (wavelength,  viewing  direction,  IOPs, 
etc)  to  be  used  in  modeling  the  Xt. 

p  oceanic  state  vector  of  retrieved  IOPs,  _1. 

p  perturbed  oceanic  state  vector  of  retrieved  IOPs,  “ 1 . 

||p||  determinate  of  p 

r  water-to-air  surface  reflectance  for  plane  irradi¬ 
ance,  dimensionless;  correlation  coefficient  when 
specifically  identified. 

R  plane  irradiance  reflectance,  in  water,  EJzyE^z), 
dimensionless. 

R[L  ratio  of  air-to-water  mean  cosines  for  downwelling 
irradiance. 


16-14 


HOGE  ET  AL  :  RTE  INVERSION  VIA  SHAPE  FACTOR 


R^l  estimate  of  ratio  of  air-to-water  mean  cosines  for 
downwelling  irradiance. 

RSR  remotely  sensed  reflectance,  in-water,  Z,M(0,  4>,  ZV 
Eocfc),  sr_1. 

RSRa  remotely  sensed  reflectance,  in  air,  Lwa(0,  $)IEoda  = 

M-da^rs*  Sr 

Rrs  remote  sensing  reflectance,  in  air,  £wa(0,  $)IEda, 
sr"1. 

RTE  radiative  transfer  equation. 

s  subscript  denoting  solar;  used  in  subscript  rs 
denotes  remote  sensing. 

5  spectral  slope  within  the  ad  model  for  CDOM+ 
detritus,  nm-1. 

t  water- to-air  radiance  transmittance,  dimensionless. 
1  water-to-air  plane  irradiance  transmittance,  dimen¬ 
sionless. 

TCA  total  constituent  absorption,  ah  m_1. 

TCB  total  constituent  backscattenng,  bbti  “l. 

v  subscript;  italicized  letter  “v”  denoting  viewing. 

V  backscattering  enhancement  factor,  defined  by 
equation  (12),  dimensionless. 

X\  shape  factors  and  related  quantities;  i  —  1,  2,  3, 
and  4,  denote  fbi  fLy  k ,  and  Ryl  respectively, 
z  depth,  m. 

OLi  model  fitting  parameters,  see  text  and  tables  for 
numerical  values. 

Oiik  model  fitting  coefficients;  a  different  set  of 
coefficients  is  needed  for  each  factor  Xt. 

(3  volume  scattering  function. 

6  Dirac  delta  function. 

A/>  uncertainty  or  perturbation  of  D 
||A/)||  determinate  of  AD. 

6*  uncertainty  or  perturbation  of  h 
1 1 6*  1 1  determinate  of  bh. 
k .(D)  ||D||  ||2>-i||,  the  condition  number  of  D 

X  wavelength,  nm. 

\b  reference  X  for  total  constituent  backscattering 
(TCB)  coefficient  model,  nm. 

\d  reference  wavelength  for  CDOM+detritus  absorp¬ 
tion  coefficient  model,  nm. 

X^  peak  wavelength  for  Gaussian  phytoplankton 
absorption  coefficient  model,  nm. 

X,  wavelength  of  observational  bands,  i  =  1,  2,  3,  . . 
nm. 

v  italicized  Greek  letter  nu  is  not  used  in  this  paper: 
see  italicized  letter  “vM  above, 
average  cosine  for  downwelling  irradiance,  in 
water,  dimensionless. 

\Lda  average  cosine  for  downwelling  irradiance,  in  air, 
dimensionless. 

<j>  azimuth  angle,  radians;  subscripts  v,  s  denote 
viewing,  solar. 

0  polar  zenith  angle  in-water  with  respect  to  +z  axis, 
radians;  subscripts  v,  s ,  and  a  denote  viewing, 
solar,  and  in-air,  respectively. 

6,  included  angle,  solar-to-viewing  direction, 
single-scattering  albedo  b/c. 
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